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THESIS  DISCLAIMER 
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available,  to  ensure  that  the  programs  are  free  of 
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without  additional  verification  is  at  the  risk  of  the  user. 
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TABLE  I    DEFINITION  OF  TERMS 

KE  =  Kinetic  Energy 

PE  =  Potential  Energy 

9  e  Large  motion  generalized  coordinate 

F  s  Large  motion  generalized  force,  applied  moment 

U  s  column  vector  for  small  motion  generalized 

coordinates, 

r  =  local  position  vector  (L  subscript  refers  to  load) 

R  s  global  position  vector  (L  subscript  refers  to  load) 

d  =  deformation  vector  of  transverse  deflections 

u  s  mass  density  of  the  link  (L  subscript  refers  to  load) 

W  =  transformation  matrix  (a  function  of  0) 

W  s  derivative  of  the  transformation  matrix  wrt  9 

W    =   residual  accelerations  term 

r 

D    =   deformation  transformation  matrix 

*    s   shape  function  matrix  derived  in  Appendix  C 

r    =   spacial  derivative  of  the  shape  function  matrix 

C     e   stiffness  matrix 

g_    =   gravitational  acceleration  vector 

I##   s   volume  integral  with  a  kernel  for  dependent  matrices 

where  ##  represents  an  index  number  from  Table  A 
K.    =   length  integral  (along  the  boom) 
H„    =   volume  integral  where  the  #  represents  an  index 
number  #  from  Table  A. 2 


TABLE  II    DEFINITIONS  OF  INTEGRAL  EXPRESSIONS 


I1;L(WTW)  =  Ju  r T  (WT  W)  rdv      (the  kernel  would  be  an 

indicated  function  of  W) 


I12(WTW)  =  Ju  rT  (WT  W)  *  dv 
I21(WTW)  =  Ju  *T  («T  W)  r  dv 
I22(WTW)  =  [u   *T  (WT  W)  *  dv 


i   -I 


T 

r  c  r  dx 


T 


Hl   "  f»  =■' 


dv 


H2   =J» 


T 
*   dv 


H4   =  K  ■£  dVl 


HK  (DTW  WDJ  ■  furrJ(D  WTWD„)  r. dv.   (given  the  tip  load  as  a 
b  K     J  L-L         K   -L   L    point  lQad^  and  the 

kernel  as  function  of 
D  and  W) 
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I.   INTRODUCTION 

A.   BACKGROUND 

The  U.S.  Navy  has  declared  it's  intention  to  enhance  its 
global,  real-time  monitoring  of  the  ocean  characteristics  and 
to  improve  its  ability  to  forecast  oceanographic  parameters 
by  means  of  a  dedicated  space  system  [Ref.  1].  The  system 
under  consideration  would  improve  the  battle  group 
commander's  ability  to  screen  his  forces  from  detection  by 
enemy  submarines,.  be  more  aware  of  severe  weather  conditions 
in  the  operating  theater  and  improve  forecasters  meteorology 
predictions . 

The  satellite  will  use  a  variety  of  special  sensors  and 
operate  from  an  820  km  polar  orbit.  A  number  of  physical 
considerations  for  this  type  of  mission  embody  the  foreground 
of  modern  structural  and  control  engineering  research. 

One  such  item  considered  necessary  to  achieve  the  Navy's 
desired  goals  is  a  Low  Frequency  Microwave  Radiometer  (LFMR). 
It  represents  a  unique  aspect  of  spacecraft  engineering  and 
control.  The  LFMR  would  provide  sea  surface  temperatures 
accurate  to  within  a  half  of  a  degree  Kelvin  with  a  surface 
resolution  capability  of  25  km.  The  LFMR  would  require  a 
large  (approximately  5  meter  diameter)  antenna  with  precise 
pointing  knowledge;  its  field  of  view  would  scan  at  either  31 
or  16  scans  per  minute.     The   width   of   scan   would   be  at 
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approximately  1400  km  centered  on  the  ground  track,  in  front 
of  or  behind  the  satellite.  The  incident  beam  would  make  a 
53.1  degree  angle  with  the  local  vertical  at  the  surface,  see 
Figure  1.  [Ref.  2:p.  3] 

The  large  ground  coverage  is  provided  by  the 
eccentrically  mounted  antenna.  A  schematic  drawing  of  the 
proposed  LFMR  sub-system  required  to  support  the  antenna  is 
presented  in  Figure  2.  The  boom  sizes  and  off-axis  angles 
vary  from  contractor  to  contractor,  those  represented  in 
Figure  2  are  not  considered  the  only  feasible  choices. 

B.   MOTIVATION 

Currently  the  Navy  .has  halted  development  on  the  Navy 
Remote  Ocean  Sensing  System  (N-ROSS),  which  was  to  be  the 
vehicle  for  the  LFMR.  Cost  considerations  and  the 
technological  development  (simulations  and  modeling  )  of  the 
LFMR  boom  assembly  and  control  scheme  have  been  suggested  as 
primary  concerns. 

This  investigation  is  motivated  by  the  limited  number  of 
techniques  currently  employed  for  the  dynamic  analysis  of 
large  flexible  space  structures.  We  will  present  the 
development  and  computer  simulation  of  an  analysis  for  a 
flexible  spacecraft  boom,  rotating  in  three  dimensions. 
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C.   LITERATURE  REVIEW 

The  approach  to  modeling  the  kinematics  and  dynamics  of 
the  rotating  flexible  space  antenna  are  not  unlike  those 
techniques  employed  in  the  field  of  robotics  necessary  for 
the  design  and  control  of  flexible  manipulators.  In  the  case 
of  flexible  manipulator  models,  the  focus  is  to  develope  the 
relationship  between  the  large,  rigid-body  motion  (large 
motion)  and  the  small  motions  due  to  the  structure's 
flexibility. 

The  most  popular  class  of  these  reduction  methods  for 
dynamic  analysis  of  structures  appears  to  employ  the  modal 
superposition  technique  [Refs.  3-7].  In  this  type  of 
analysis  a  determination  must  be  made  to  limit  the  number  of 
contributing  mode  shapes  and  then  combine  those  that  will 
contribute  in  a  way  which  yields  a  system  of  equations  of 
motion. 

Book  [Refs.  4,5]  utilizes  this  modal  technique  to  model  a 
system  with  flexible  kinematics.  He  utilizes  recursive 
Lagrangian  kinematics  and  large  and  small  motion  coordinates. 
The  resulting  equations  of  motion  are  an  accurate  model 
describing  large  and  small  motion  and  the  coupling  between 
them.  However  the  result,  a  set  of  coupled,  non-linear 
second  order,  ordinary  differential  equations,  proved 
computationally  intensive,  and  therefore  expensive  to  solve. 

Sunada  and  Dubowsky  [Ref.  8]  on  the  other  hand  expressed 
the  kinematics  and  dynamics   of  the   manipulator  in   terms  of 
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"small  perturbations  about  some  known  nominal  motion  of  the 
system".  The  effect  of  the  coupling  between  the  large  and 
small  motion  is  incomplete  when  the  large  motion  is  assumed 
nominal  and  the  small  motion  superimposed  upon  it. 

Chang's  development  of  the  Equivalent  Rigid  Link  System 
ERLS  [Ref.  9],  describes  the  large  motion  kinematics  of  the 
system  by  the  ERLS  and  the  small  motion  deformations  relative 
to  it.  The  application  of  finite  element  techniques  and 
Lagrangian  dynamics,  produces  two  sets  of  coupled,  non- 
linear, ordinary  differential  equations  of  motion.  The 
application  of  the  ERLS  formulates  these  sets  of  equations; 
one  set  for  large  motions  and  one  set  for  small  motions.  The 
set  of  large  motion  equations  are  non-linear  in  both  the 
large  and  small  motion  variables.  The  set  of  small  motion 
equations  are  linear  in  the  small  motion  variable  and  non- 
linear in  the  large  motion  variables. 

The  characteristics  of  these  equation  sets  permit  a 
solution  by  numerical  integration  with  a  variety  of 
techniques  possible.  The  integration  technique  developed  by 
Chang,  the  Sequential  Integration  Method,  is  straightforward 
and  relatively  easy  to  implement.  Chang's  model  offered  a 
complete  representation  of  the  dynamics  of  flexible 
manipulators  as  well  as  one  that  could  be  efficiently  solved 
using  the  Sequential  Integration  Method. 

Simo  and  Vu-Quoc  [Ref.  10]  present  an  approach  which 
appears  well   suited  for  design  but  the  choice  of  generalized 

15 


coordinates  may  require  a  complicated  observer  to  relate  the 
computational  prediction  of  the  generalized  coordinates  with 
observable  reference  quantities  of  the  controller 
environment. 

Previous  work  at  the  Naval  Postgraduate  School  utilizing 
the  ERLS  technique  validated  the  procedure  in  a  planar 
application  [Refs.  11,12].  Both  Petroka  and  Gannon  applied 
Chang's  ERLS  utilizing  separate  mode  shape  response  matrices. 
The  agreement  between  the  assumed  modal  responses  and  the 
experimentally  observed  response  justifies  this  first  effort 
at  applying  the  ERLS  to  model  the  dynamic  response  of 
flexible  booms  in  the  spacecraft  applications. 

D.   OUTLINE 

The  nonlinear  dynamic  response  of  a  spacecraft  boom 
structure  member  is  predicted  by  the  numerical  integration  of 
the  equations  of  motion.  In  this  work,  the  spatial  finite 
element  discretization  of  the  boom  and  the  application  of  an 
assumed  polynomial  mode  shape  provides  an  approximate 
solution. 

The  modeling  process  has  been  applied  through  the  ERLS 
utilizing  Lagrangian  Dynamics  to  obtain  the  system  equations 
of  motion  for  numerical  integration. 

Chapter  Two  of  this  work  reviews  Chang's  development  of 
the  ERLS  and  in  Chapter  Three  the  ERLS  is  applied  to  a 
rotating  flexible   boom  and   a  variety   of  simulation  results 

16 


are  presented.  The  computer  simulation  code  and  methodology 
are  introduced  in  Chapter  Four.  The  results  of  this  work 
illustrate  both  the  need  for  controlling  the  flexibility  and 
the  versatility  of  the  ERLS  solution  technique.  A  review  of 
this  work  and  suggested  areas  for  development  is  presented  in 
Chapter  Five. 
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II.   DYNAMIC  MODEL  FORMULATION  FOR  ft  SPACECRAFT 
SINGLE-ELEMENT  FLEXIBLE  BOOM 


A.   KINEMATIC  RELATIONSHIPS  OF  THE  EQUIVALENT  RIGID  LINK 
SYSTEM 

The  introduction  and  development  of  the  Equivalent  Rigid 
Link  System,  separates  the  kinematics  of  the  flexible  body 
into  large  and  small  motions.  The  large  motions  are  modeled 
with  the  characteristics  of  rigid  bodies  while  the  small 
motions  are  modeled  by  finite  element  techniques  and  defined 
relative  to  the  position  of  the  equivalent  rigid  bodies. 

The  small  motions  resulting  from  loading  flexure  are 
measured  in  the  ERLS  relative  to  the  local  coordinate  system 
for  each  equivalent  rigid  body.  The  absolute  position  of  a 
point  on  a  deformed  body  is  obtained  by  coordinate 
transformations  which  utilize  the  characteristics  of  joint 
variables  as  in  the  methods  of  Hartenberg  and  Denavit  [Ref. 
13].  Lagrangian  dynamics  are  used  to  obtain  the  equations  of 
motion  and  therefore  time  derivatives  of  the  absolute 
position  will  be  necessary  in  the  development  of  the 
Lagrangian  expression. 

The  ERLS  model  accounts  for  the  change  in  location 
experienced  by  each  point  along  a  flexible  member.  Discrete 
points  or  nodes  are  chosen  and  their  relative  position 
determined  while  locations  of  intermediate  points  are 
assumed  according   to  a  mathematical  function.   The  choice  of 
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nodes  and  mathematical  function  determines  the  theoretical 
shape  of  the  flexed  member. 

The  Finite  Element  Method  (FEM)  is  used  to  descretize  the 
flexible  member  assigning  the  nodes  and  enforcing 
compatibility  between  elements.  The  choice  of  shape  function 
and  node  position  must  be  consistent  with  the  accuracy 
desired.  The  ERLS  will  accommodate  closed  form  natural  mode 
shape  functions  or  any  order  of  polynomial.  In  this  work  for 
example  a  cubic  polynomial  shape  function  was  chosen  to 
represent  the  bending  displacement  of  the  flexible  boom.  The 
nodal  positions  chosen  here  are  the  end  points  of  the 
flexible  member.  The  selection  of  a  polynomial  shape 
function  and  end  point  nodal  positions  are  consistent  with 
previous  works  which  have  experimentally  validated  their 
selection.  [Refs.  9,11,123. 

The  kinematics  of  the  ERLS  separates  the  global  motions 
into  the  sums  of  large  and  small  motions.  The  large  motion 
is  described  by  the  equivalent  rigid  link  and  the  small 
motion  due  to  deformations  of  the  link  are  described  relative 
to  the  rigid  link. 

The  characteristics  of  rigid  body  dynamics  is  essentially 
classical  dynamics.  In  this  application  the  Hartenberg- 
Denavit  4x4  transformation  matrix  is  used  to  completely 
describe  the  positions  and  velocities  of  each  point  in  the 
system.  Appendix   A  develops  these  kinematic  relationships  in 
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a  brief  presentation,  a  complete  review  of  the  kinematics  of 
the  ERLS  is  available  [Ref  9:pp. 17-50]. 

B.   KINETIC  RELATIONSHIPS  OF  THE  EQUIVALENT  RIGID  LINK 
SYSTEM 

The  formulation  of  the  Lagrange  equations  to  derive  the 
equations  of  motions  is  straightforward  and  systematic.  This 
familiar  approach  is  helpful  in  the  complex  analysis  of  this 
type. 

The  selection  of  generalized  coordinate  systems  is 
accomplished  indirectly  by  the  separation  of  the  global 
motion  into  large  and  small  motions.  Formalizing  this  choice 
assigns  the  large  motion  generalized  coordinate  as  theta  (9), 
which  represents  a  rotation  and  the  small  motion  generalized 
coordinates  (U)  a  column  vector  composed  o#f  the  translations 
and  rotations  at  the  boom  tip. 

The  translations  from  the  nodal  position  on  the  ERLS  to 
the  deformed  position  on  the  actual  flexible  member  (x,y,z) 
will  become  the  first  three  small  motion  generalized 
coordinates.  The  slopes  at  the  node  or  the  rotational  angles 
of  flexible  member  about  the  axis  perpendicular  to  the  planes 
of  displacement  of  the  flexible  member  (U')  are  included  as 
the  last  three  small  motion  generalized  coordinates. 

The  transformation  from  the  ERLS  to  the  deformed  boom 
position  is  accomplished  by  the  mode  shape  matrix.  The 
derivation  of  the  mode  shape  matrix   is  briefly   presented  in 
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Appendix  B.  The  choice  of  nodal  position  and  mode  function 
determines  the  accuracy  of  the  simulation  response. 

In  this  representation  a  cubic  polynomial  function  was 
chosen  as  the  shape  function  for  bending  and  a  linear 
function  for  the  extension  and  torsion.  This  choice  provides 
the  mathematical  representation  of  our  first  three  small 
motion  generalized  coordinates,  and  simplifies  the  required 
mathematics  to  obtain  the  functional  representation  of  the 
remaining  three  small  motion  generalized  coordinates. 

The  system  kinetic  energy  is  the  sum  of  the  kinetic 
energies  of  each  member  of  the  system,  this  would  include  the 
flexible  boom  and  the  tip  load  (if  any).  The  system 
potential  energy  results  from  the  elastic  strain  energy  of 
the  flexible  member,  and  the  application  of  gravitational 
body  force  (if  any).  Generalized  forces  represent  the 
application  of  applied  forces  and  damping  forces.  Appendix  C 
includes  the  development  of  the  equations  of  motion  from  the 
Lagrange  Equations  for  the  case  of  end  point  nodal  positions 
and  a  polynomial  shape  function  for  six  small  motion 
generalized  coordinates. 

The  manipulations  of  the  Lagrange  equations  to  obtain  the 
equations  of  motion  entail  extensive  effort,  mathematical 
dexterity,  and  fierce  determination  but  the  reward  is  a  set 
of  two  non-linear,  coupled,  second-order,  ordinary 
differential  equations.    Though   both  sets   are  coupled,  one 
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describes  the  large  motions  and  the  other  describes  the  small 
motions,  they  are  presented  below. 

M    e'  +  M   U-  =  F  (2-1) 

qq        qn       q 

M     6*    +  M    U%  G   U  +  K   U  =  f       (2-2) 
nq        nn      n      n      n 
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lxl  effective  inertia  matrix  for  large  motions 

1x6  coupled  effective  inertia  matrix  of  the 

small  motion  on  the  large  motion 

lxl  load  vector  for  the  large  motion 

6x1  coupled  effective  inertia  matrix  of  the 

large  motion  on  the  small  motion 

6x6   effective    inertia   matrix    for    small    motions 

6x6  gyroscopic  matrix 

6x6  stiffness  matrix 

6x1  load  vector  for  the  small  motion 

generalized  coordinate  of  the  large  motion 

6x1  generalized  coordinate  vector  for  small 

motion  deformations 


C.   NUMERICAL  INTEGRATION  OF  THE  EQUATIONS  OF  MOTIONS 

The  application  of  numerous  numerical  integration  schemes 
is  possible  in  the  solution  of  the  equations  of  motion.  The 
primary  consideration  in  the  selection  of  the  integration 
technique  is  the  determination  of  the  time  step  necessary  to 
integrate   the   equations   of   motion   and   obtain   numerical 
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stability.  In  this  work,  several  integrations  techniques  were 
tested  i.e.:  Runga-Kuta  5th  order  (both  fixed  and  variable 
time  step),  Adams -Moul ton  (variable  time  step),  and  Gear. 
While  all  of  these  techniques  achieved  numerical  stability 
the  computational  time  required  to  obtain  adequate  simulation 
time  was  excessive. 

The  formulation  of  the  equations  of  motion  in  a  manner 
similar  to  equations  (C-36)  and  (C-37),  permits  the 
application  of  the  Sequential  Integration  Method  (Ref.  9:pp 
79-911.  In  this  technique  an  implicit  integration  scheme  is 
employed  to  integrate  the  small  motion  equation,  and  an 
explicit  technique  is  used  to  integrate  the  large  motion 
equation. 

The  implicit  technique  guarantees  numerical  stability  in 
the  case  of  linear  equations  and  hence  an  iterative  technique 
can  be  avoided.  The  results  of  the  small  motion  integration 
can  be  provided  to  the  large  motion  equation  which  can  then 
be  integrated  explicitly  in  a  predictor-corrector  fashion. 
The  choice  of  time  step  in  the  case  of  the  explicit 
integration  can  be  governed  by  solution  accuracy  vice  small 
motion  frequency  response.  The  decided  advantage  to  this 
approach  is  the  computational  time  savings. 
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III.   APPLICATION  OF  THE  ERLS  TO  A  SPACECRAFT 
SINGLE  ELEMENT  FLEXIBLE  BOOM 


A.   THE  LOW  FREQUENCY  MICROWAVE  RADIOMETER  (LFMR) 

The  low  frequency  microwave  radiometer  (LFMR)  which  would 
be  based  on  a  space  platform  of  the  N-ROSS  type  is  composed 
of  a  reflector,  an  upper  reflector  boom,  a  lower  reflector 
boom,  and  a  spacecraft  boom.  These  booms  would  all  be  hinged 
together  and  owing  to  the  requirements  of  space  systems, 
would  be  developed  from  light  weight  materials.  The  dynamic 
modeling  of  such  a  structure  has  been  the  primary  focus  of 
the  work  presented  here,  other  applications  of  the  ERLS  in 
space  technology  do  however  exist. 

The  orientations  of  the  microwave  receiving  dish  on  board 
requires'  precise  pointing  accuracy  which  must  be  well 
controlled  for  the  transmission/reception  of  signals.  The 
accuracy  of  the  pointing  mechanism  is  subject  to  the  position 
error  determinations  which  in  turn  are  vulnerable  to  the 
accelerations  of  orbital  adjustment,  antenna  spinning  and  the 
inherent  light  weight  materials  of  construction.  For  these 
reasons  the  control  of  vibrations  is  essential. 

The  vibrational  response  or  deformation  resulting  from 
the  application  of  these  forces  must  be  observable  to  the 
controller.  The  choice  of  the  small  motion  coordinates 
previously  outlined  would  make  these  deformations  available 
during   start   up,   regular   operation,   and   orbital  station 
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keeping  adjustments.  Since  these  values  are  directly 
measurable  this  choice  could  simplify  the  controller  observer 
requirements  while  ensuring  the  mission  specifications  are 
met. 

B.   SIMULATION  OBJECTIVES 

In  this  work  the  multi-bodied  boom  has  been  simplified  to 
a  single  member.  The  deformations  resulting  from  the 
structural  flexibility,  have  been  assumed  small  and  hence 
small  angle  assumptions  have  been  implemented  in  the  model. 
The  application  of  body  forces  and  internal  forces  have  been 
limited  to  gravity,  and  a  driving  torque.  No  damping  effects 
have  been  implemented,  nor  contributions  to  deformations  from 
thermoelastic  effects  upon  entering  or  leaving  the  eclipse. 
The  incorporation  of  these  concerns  is  possible  in  a  similar 
model  but  beyond  the  scope  of  the  work  presented  here. 

The  exact  material  of  the  boom  would  likely  be  a  graphite 
epoxy  composite  with  specific  material  properties,  however 
the  simulation  code  utilizes  the  material  considerations  of 
pure  aluminum  to  permit  the  intuitive  interpretations 
associated  with  a  more  familiar  structure  material.  Specific 
material  considerations  can  be  easily  changed. 

The  ultimate  goal  of  this  work  was  to  complete  the  first 
stages  of  the  simulation  modeling.  To  develope  the  solution 
techniques   and   the   basic   simulation  code  (Appendices  A-D) 
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which  can  later  be  expanded  to  eventually   include  multi- 
elements  in  the  satellite  model. 

C.   DESCRIPTION  OF  THE  COMPUTER  SIMULATION  CODE 

The  IBM  simulation  language  DSL  (Dynamic  Simulation 
Language)  was  originally  chosen  to  model  the  equations  of 
motions.  The  variety  of  integration  schemes  and  FORTRAN 
compatibility  were  the  over  riding  reasons  for  this  choice. 
Gannon's  work  in  two  dimensions  indicated  a  preferred 
integration  technique  for  the'  Adams-Moulton  method  [Ref. 
12:p.  25],  additionally  DSL  proved  capable  of  simulating 
Gannon's  system  model. 

While  integration  methods  were  not  specifically  the  focus 
of  this  work,  as  indicated  previously,  the  choice  of 
integration  technique  became  dependant  upon  the  time  step 
requirements  for  numerical  stability. 

The  simulation  code  was  developed  with  modular  FORTRAN 
subroutines  then  attached  to  a  short  initialization  program. 
From  the  three  dimensional  development  in  Appendix  C  the 
Lagrange  equations  can  be  formulated  as  equations  (2-1)  and 
(2-2)  and  then  expanded  in  the  form  of  equations  (C-36)  and 
(C-37).  Where  the  formulation  in  this  manner  lends  itself  to 
the  Sequential  Integration  Method. 
The  resulting  column  vector: 

[9uuuuuu]  (3-1) 

x        y      z      xx      yy      zz 


26 


represents  the  large  motion  acceleration  8  and  the  rate  of 
local  deformation  in  position  and  slope.  Where  the  x  y  and  z 
subscripts  represent  the  local  coordinates  and  the  double  x  y 
and  z  subscripts  represent  the  rotations  about  these  axis. 

Each  subroutine  of  the  simulation  code  was  developed  to 
specifically  determine  a  portion  of  the  single  element  three 
dimensional  analysis  and  due  to  the  nature  of  modular 
programing  the  expansion  of  this  code  to  a  multi-body  problem 
can  be  accommodated.  The  logic  flow  of  the  coding 
(Figure  15)  is  provided  to  assist  in  the  future  development 
of  multi-body  simulations. 

The  simulation  code  can  be  viewed  from  three  levels; 
LEVEL  I  (an  overview)  separates  the  code  into  the  primary 
focal  points  of  INITIALIZATION,  FORMULATION,  INTEGRATION, 
REFORMULATION,  and,  OUTPUT.  The  reformulation  is  a  technique 
to  incorporate  the  inertia  coupling  between  the  large  and 
small  motion  equations. 

LEVEL   II   can   be   viewed   as  the  FORTRAN  formulation  of 
coefficients  matrices,  and  sequential  integration  to  solve 
for  9,  6,  9*,  U,  U,  and  5". 

LEVEL  III  is  a  listing  of  the  FORTRAN  utilities  sub- 
programs required  for  the  matrix  formulation  and  manipulation 
in  LEVEL  II.  Each  of  the  units  in  LEVEL  II  may  access  the 
routines  in  LEVEL  III. 
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D.   SIMULATION  RESULTS 

The  computer  simulation  code  was  run  with  a  variety  of 
initial  conditions  these  conditions  were  chosen  to  present 
the  aspects  of  the  boom  flexibility  and  the  impact  this 
flexibility  would  have  on  the  requirements  for  end  point 
control.  The  simulation  output  has  been  presented  in  local 
displacement  values,  where  these  relative  values  represent 
the  conditions  at  the  boom  tip,  the  point  mass  location. 

The  simulation  results  presented  in  Appendix  E  are 
considered  typical  of  results  available  from  the  simulation 
code.  They  represent  the  single  large  motion  generalized 
coordinate  and  each  of  the  six  small  motion  generalized 
coordinates.  The  intent  of  these  figures  is  to  provide  the  a 
physical  feeling  for  the  plant  behavior. 

The  results,  Figures  3-7,  reflect  the  slowing  effects  on 
the  large  motion  generalized  coordinates  from  the  point  mass. 
The  small  motion  effects,  Figures  8-14,  and  the  coupling 
between  the  large  and  small  motion  coordinate  systems,  is 
also  clearly  observable.  Additionally  the  results  indicate 
the  effects  of  boom  elevation  angle  and  corresponding 
responses  at  four  elevations,  zero  through  60°. 
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IV.   SUMMARY 

A.   CONCLUSIONS 

This  investigation  is  motivated  by  the  needs  of  the 
United  States  Navy  to  be  able  to  accurately  predict  the 
Meterological  phenomenon  within  an  operations  theater.  To 
this  end  the  techniques  currently  employed  for  the  dynamic 
analysis  of  large  flexible  space  structures  must  be  improved 
and  expanded.  We  have  presented  the  development  and  computer 
simulation  of  an  analysis  for  a  flexible  spacecraft  boom, 
rotating  in  three  dimensions  with  flexible  modes  in  bending, 
extension  and  torsion. 

The  nonlinear  dynamic  responses  of  the  structure  member 
are  predicted  by  the  numerical  integration  of  the  equations 
of  motion.  In  this  work,  the  spatial  finite  element 
discretization  of  the  boom  structure  and  the  application  of 
an  assumed  third  order  modal  response  provide  an  approximate 
solution  to  the  equations  of  motion,  one  which  lends  itself 
to  the  problems  associated  with  controller  requirements. 

The  modeling  process  has  been  applied  through  the 
Equivalent  Rigid  Link  System,  utilizing  Lagrangian  Dynamics 
to  obtain  the  system  equations  of  motion  for  numerical 
integration.  The  effectiveness  of  the  integration  schemes 
available  for  numerical  integration  of  this  type  of  problem 
have  been  briefly  presented. 
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The  justification  for  the  ultimate  selection  of  a 
Sequential  Integration  Method  was  discussed.  Areas  of  the 
computer  code  which  must  be  developed  in  order  to  extend  the 
model  to  multiple  bodies  and  eventually  to  the  controller 
design  and  application  were  also  presented. 

B.   RECOMMENDATIONS 

Areas  which  remain  to  be  investigated  utilizing  the 
simulation  code  include: 

1.  The  significance  of  torsion  and  extension  in  the 
total  system. 

2.  The  significance  of  the  way  in  which  torque  is 
applied  to  the  boom,  for  example;  impulse,  ramp  or  step 
wise  changes. 

3.  The  significance  of  the  material  properties  of  the 
model . 

4.  The  significance  of  the  rotational  velocity. 

The  single  largest  concept  remaining  to  be  addressed  by 
this  boom  model  is  the  additional  members  to  be  joined  to  the 
driven  boom.  The  addition  of  these  members  requires 
additional  application  of  the  ERLS  modeling  scheme. 
Implementation  requires  compatibility  for  both  members  at 
the  connecting  joints. 

The  long-term  goal  of  this  research  is  the  design  and 
control  of  a  large  flexible  antenna-boom  assembly  for  a  space 
based   meterological   data   acquisition   system.    Simulation 
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studies  need   to  continue   and  their  scope  should  be  expanded 
to  multi-body  analysis. 

The  control  system  requirements  must  be  investigated  and 
incorporated  into  the  antenna  system  model.  Finally  a 
physical  model  must  be  developed  to  validate  the  simulation 
results . 
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APPENDIX  A 

DERIVATION  OF  THE  KINEMATICS 
REQUIRED  IN  THE  EQUATIONS  OF  MOTION 


In  Appendix  C  we  presented  the  derivation  of  the 
equations  of  motion  for  the  single-element,  flexible  boom. 
In  the  derivation  of  these  equations  we  represent  the 
position  of  any  point  in  the  deformed  member  globally  by 
applying  either  a  coordinate  transformation  matrix  W  or  a 
deformation  transformation  matrix  D  and  in  the  determination 
of  the  small  motions  the  application  of  both  transformations 
is  required. 

In  order  to  more  accurately  represent  the  coordinate 
transformation  matrix  W,  we  will  present  the  Hartenberg- 
Denavit  transformation  matrix  A,  and  demonstrate  the 
technique  for  obtaining  the  matrix  W  and  it's  derivatives. 
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(A-l) 
The   theory   and   derivation   of   the   Hartenberg-Denavit 
matrix  formulations  are  defined  in  several   of  the  references 
cited  in  Chapter  One  this  work  [Refs.  9,11-13].   The  intent 
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here   is   to   present   the   application  of   this  technique. 
Variables  variables  used  in  the  H-D  transformation  matrix  are 
called  joint  variables  and  these  are  defined  below. 

9  =  the  rotation  of  the  link  about  the  local  z  axis 

a  =  the  length  of  the  flexible  link 

d  =   the  z  axis  displacement  between  consecutive 

coordinate  systems 
a  =  the  angular  rotation  between  consecutive  local  z 

axis 

In  the  application  of  the  H-D  matrix,  the  rotation  angle 
6,  of  the  member,  is  always  about  the  local  z  axis.  To 
obtain  this  orientation  we  have  defined  two  sets  of  local 
coordinates  which  will  require  two  coordinate 
transformations.  The  first  transformation  is  a  rotational 
one,  we  transform  the  global  z  axis  to  the  local  y  axis, 
(this  is  a  rotation  of  a  radian  about  the  x  axis).  The 
transformation  matrix  Al  represented  in  equation  (A-2), 
provides  this  rotational  transformation. 
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The  second  local  coordinate  system  is  defined  at  the  end 
of  the   link,  and   therefore  we   have  included  the  effects  of 
translation.    We  have  represented  this  as  A  .  it  is  defined 
by  equation  (A-3). 
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Having  defined  the  two  H-D  matrices  required  in  the 
simulation,  we  can  represent  the  coordinate  transformation 
matrix  W. 


W  =  I   A2 


(A-4) 


The  definition  of  W  lends  itself  toward  the  application 
of  the  chain  rule  in  the  determination  of  the  derivatives  of 
W.  By  the  choice  of  the  H-D  transformation  matrix  we  can 
easily  obtain  the  derivatives  of  each  A  matrix  as  follows: 


dq 


A  =  Q  A 


(A-5) 


Where  the  derivative  is  taken  with  respect  to  the 
generalized  coordinate,  in  this  application  q  would  be  the 
angle  theta  and  Q  is  defined  by  equation  (A-6), 
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Q  = 
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(A-6) 


If  we   applied   the   unique  conditions  of  the  spacecraft 
rotating  boom  we  can  set  the  change  in  q2  as  zero  and 

eliminate  the  concern   for   order   in   multiplication   of  the 
three  matrices  which  compose  the  derivative  of  the  W  matrix. 


and: 
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For  the  time  derivatives  we  obtain  (with  6   constant): 


W  -   ^-  = 


dW   _9_W    _96_   _9W    • 
dt  "  ge^      at  "  d^i        1 


(A-9) 


•'  =  ar   *  ■  ar  QAi*iA2  =  QAi5iliA2  +  QAiK  A2 


(A-10) 


In  expression   (A-10)  we   can  consolidate   the  product  of 

each  group  of  terms  on  the  extreme  right  hand  side  and 

•  2 
rename  them  the  residual  acceleration  W  9.  and  the  large 

r   l 

■  ■ 
motion  acceleration  We  .   In  following  derivations 

subscripts  for  0  have   been  omitted   since  there   is  only  one 

large   motion   generalized   coordinate  capable  of  motion,  the 

derivatives  of  the  other  would  be  zero. 
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Defining  the  deformation  transformation  matrix  as 
equation  (A-ll)  will  accommodate  three  dimensional 
deformation. 
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(A-ll) 


The  partial  derivative  of  D  with  respect  to  each  of  the 
small  motion  generalized  coordinates  results  in  the  same 
matrix  as  the  partial  of  the  time  derivative  of  D  with 
respect  to  the  time  derivative  of  the  small  motion 
generalized  coordinates,  equation  (A-12).  In  the  following 
development  the  subscript  K  refers  to  the  particular  small 
motion  generalized  coordinate  of  interest. 
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APPENDIX  B 
DERIVATION  OF  THE  MODE  SHAPE  RESPONSE  MATRIX  * 

The  derivation  of  the  mode  shape  matrix  has  been 
presented  previously  [Refs.  9,11,12].  The  derivation 
presented  here  is  abbreviated  and  included  only  as  a  brief 
description  of  the  techniques  involved,  for  complete 
derivation  the  reader  is  encouraged  to  review  the  previous 
works . 

Representing  the  relative  position  changes  due  to 
deformations  by  d_  we  have  assumed  small  deformations,  no 
warping  after  torsion  and  all  deformations  u,  v,  w,  and  0  are 
only  functions  of  x  (along  the  boom),  the  yz  plane  remains 
plane  after  torsion  and  deformations  only  vary  along  the  x 
direction. 
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(B-l) 


Applying  simple  beam  theory,  we  have  assumed  that  the 
rotary  effect  and  the  lateral  shear  deflections  are 
negligible. 
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We  have  chosen  a  polynomial  form  to  represent  the  shape 
function  for  the  displacement  of  the  single  element  flexible 
boom  as  shown  below. 


ul  =  a0  +  alx 

U2  =  b0  +  blX  *  b2x2  +  b3x3 

U3  =  C0  +  C1X  +  C2x2  +  C3x3 


u„  =  a 


uc  =  b,  +  2b0x  +  3b0x' 
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u,  =  c,  +  2c_x  +  3c„x' 
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We  can  now  apply  the  boundary  conditions  of  zero 
displacement  and  slope  at  the  base  where  x  is  equal  to  minus 
the  boom  length,  L,  and  the  sign  is  consistent  with  the 
convention  for  the  local  coordinate  system. 

Substituting  the  boundary  conditions  into  the  shape 
functions  gives: 
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u. 


u, 


u. 


u, 


u, 


u, 


u. 


u, 


u, 


u 


u, 


u, 


-L) 
-L) 

-L) 

-L) 

-L)  = 

-L)  - 


aQ  +  a.x  =  0 

bQ  +  bxx  +  b2x2  ♦  b3x3  =  0 

2       3 
cQ  +  c.x  +  c«x   +  c3x   =  0 


=  a. 


b.  +  2  b«x  +  3  b3x 


c-  +  2  c«x  +  3  c3x 


0) 
0) 
0) 
0) 
0) 
0) 


=  a. 


=  b. 


=  c, 


=  a. 


=  b. 


=  c. 


(B-9) 

(B-10) 

(B-ll) 

(B-12) 
(B-13) 

(B-14) 

(B-15) 
(B-16) 
(B-17) 
(B-18) 
(B-19) 
(B-20) 


The  coefficients  a«,  a3,  b«,  b3,  c~,    c3,  can  now  be 
solved  for  by  substitution  and  arranged   in  the   shape  matrix 
representation  as; 


*   = 


"21  " 
0   a 


32 
0    a 


34 


36 


43  a44  a45  ° 


(B-21) 


where:  R  =  — =-—   x,  the  variable  of  integration,  L  the  length 

over  which  the  integration  will  take  place  (the  length  of  the 

flexible  member). 
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21 


=   1  +  R 


32 


l43 


34 


l44 


l45 


36 


3      2 
-  2R   -  3R   +  1 


3      2 
-  2RJ  -  3R^  +  1 


-  Z  (1  +  R) 


-  Y  (1  +  R) 


LR3  +  2LR2  +  LR 


LR3  +  2LR2  +  LR 


For  the  potential  energy  it  will  be  necessary  to  develop 
the  expression  for  u',  \i'l,    \i'l,    u',u4,  u',  accomplished  by 
differentiating  the   shape  function  equations  (B-9)  -  (B-20) 
and  then  solving  the  simultaneous  equations  a  second   time  to 
produce  the  modified  shape  function  matrix  which  is  presented 
below: 


r   = 


cll  ° 


0 
c, 


'32 
0    0    0 


~23 
0    0 


25 
0    c 


36 
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(B-22) 
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11 


1 

L 


23 


12R 

r  2 


25 


6R 


32 


23 


'36 


25 


C44    =      Cll 
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APPENDIX  C 

DERIVATION  OF  THE  EQUATIONS  OF  MOTION  FOR  THE 
SINGLE-ELEMENT,  FLEXIBLE  BOOM 


In  derivations  of  this  type  additional  definitions  are 
made  to  simplify  the  terms  involved  in  the  Lagrange 
equations.  A  table  of  terms,  for  both  variable  names  and 
integral  definitions  (Tables  I,  and  II  respectively)  have 
been  included  for  this  purpose. 

In  the  case  of  a  rotating  flexible  boom  the  large  and 
small  motions  can  be  separated  by  application  of  the  ERLS. 
The  large  motion  generalized  coordinate  is  defined  as  8  and 
the  small  motion  generalized  coordinates  as  U.  In  this 
application  U  represents  a  column  vector  of  six  entries.  We 
can  imagine  this  vector  of  generalized  coordinates  as 
composed  of  two  sets.  The  first  set  ul,  u2,  and  u3  represents 
the  translational  displacement  of  the  end  point  in  the  local 
x  y  and  z  directions  when  deformed.  The  second  set  u4,  u5, 
and  u6  represents  the  rotation  or  slope  of  the  flexed  member 
at  the  deformed  end.  In  accordance  with  this  description  the 
vector  of  generalized  coordinates  can  be  defined  by  equation 
(C-l) . 

T 


"  e    =  [9   Ul   U2   U3   U4   U5   U6 


(C-l) 
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These  generalized  coordinates   will   be   used    in   the 
application  of    the   Lagrange   equations   to   develope   the 
equations  of  motion  for  both  the   large  and   small  coordinate 
systems.    The  Lagrange  equations  are  defined  in  equations 
(C-2)  and  (C-3). 

d  r  3KE]    3KE      3PE      „  {r   0. 


d    TaKEl    _    3KE      +      3PE      =  _3 

L  3fiJ 


ar  "37      "W   +   ttj 


Small  motion  generalized  forces  are  zero  in  this  example, 
this  is  consistent  for  a  system  without  damping.  Utilizing 
the  ERLS  we  can  define  the  global  position  of  the,  end  point 
in  terms  of  the  local  position  vector  (r)  and  a  deformation 
vector  (d)  equation  (C-4).  Similarly  the  global  position  of 
the  tip  load  is  obtained  by  transforming  a  local  position 
vector  for  the  load,  vector  (rL),  first  by  a  deformation 
transformation  matrix  D  and  then  by  a  coordinate 
transformation  matrix  W. 

R  -  W  (r  +  d)  (C-4) 

El=  WDr^  (C-5) 

In  order  to  completely  represent  the  coordinate 
transformation  matrix  we  have  presented  the  Hartenberg- 
Denavit    transformation   matrix   A,   and   demonstrate   the 
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technique  for  obtaining  the  matrix  W  and  it's  derivatives  in 
Appendix  A.  Thus  we  continue  with  the  development  of  the 
equations  of  motions,  not  encumbering  the  derivation  of  these 
equations  with  the  derivations  of  others.  Also  presented  in 
Appendix  A  is  the  definition  of  the  deformation 
transformation  matrix  D. 

Continuing;  we  present  the  derivative  with  respect  to 
time  for  the  global  positions  in  order  to  obtain  the  kinetic 
energy  expressions. 

R  =  W(r  +  d)  +  Wd  (C-6) 

L  =  WDr, +  WOr,  (C-7) 


The  kinetic  and  potential  energy  expressions   equations 

(C-8)   and   (C-9)   follow.  In   this   development   we   have 

separated  the  large  and  small  motion  energy  contributions  and 
will  present  them  separately. 

For  the  large  motion: 


KE 


=  \    [  Ju  RT  Rdv  +   fu  §*  ^  dvj  (C-8) 


link  load 

volume         volume 


PE 


=  \   JuTrTC  rU  dx  -   u  B.T  g  dv  -   uL  r£  g  dvL       (C-9) 


link  link  load 

length  volume         volume 
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Now  redefining  d  in  terms  of  the  shape  matrix  *  (again 
this  derivation  is  presented  in  Appendix  B  to  preserve  the 
continuity  here).    We  can  rewrite  equations  (C-9)  as: 

KE  =  fu[i[(r+*U)TWTW(r+*U)  +  (*U)TWTW*U]  +  (r+*U) TWTW*U]  dv 

+  ^[^.[(Di^rtrWtDr^  +  (Dr^lTWtDr^  ]  +  (Dr^&TWCDr^  ]  dvL 

(C-10) 

Continuing   in   the   development   we   first   express   the 

derivative   of   the  kinetic   energy   with  respect  to  the  time 

■ 
rate  of  change  of  the  large   motion  coordinate,   0,  equation 

(C-ll),   and   then   obtain   the   time   rate  of  change  of  this 

expression,  equation  (C-12). 

£2£  =  fw[(r+*U)Tw£w(r+*U)  +  ( r+*U)TWQW*5]  dv  + 

JuL[(Dr£Tw£w(Dr£   +    (Dr£TwJw( Dr^)    ]    dvL  (C-ll) 

^[-^f-]    =  Ju[(*U)TwJw(r+*U)    +    (r+*U)T(WQW+WQW')(r+*U)    + 

8 

(r+*U)T(WQW)(*U)    +    (*G)T(WqW)(*U)    + 

(r+*U)T(wJw+WQW)(*U)    +    (r+*U)T(WQW)(*u')]    dv    + 

JuL  [(6r^TwJw(Dr^   +    (Dr^W^M   +   wJw'MDr^   + 

(Dr^TW^W(5r^   +    (fir^W^WtSr^   +    (Dr^W^tB'r^    + 
(Dr£T(wJw    +    wJwMDr^]    dvL  (C-12) 
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Finally  we  express  the   derivative  of   the  kinetic  energy 
with  respect  to  the  large  motion  position  as  0. 

^|  =Ju[j[(r+*U)T(WQW  +  WTWe)(r+*U)  +  (*0)T(WqW  +  WT«e)(*U)] 

+    (r+*U)T(WQW+flTWe)(*U)]    dv    +   Ju[|[(Dr£T(flQW    +   flTWQ)(Dr£   + 

(Sr^fWgW    +   WTWe)(Dr^]    f    (Dr^Tw^W    +   WTWQ)(Dr£]    dvL 

(C-13) 


d4C"^^    "  t!  =  J»[^+*u)T^wJw')(r+*U)    +    2(r+*U)TWQW(*U).    + 

(r+*U)TWQW(*u")]  dv  f  JVL[(Dr£TWQW*(Dr£  + 
2(Dr^TWQW(Dr^  +  ( Dr£TwJw(S'r£]  dvL       (C-14) 


The  small  motion  equation  is  handled  in  a  similar 
fashion,  however  it  is  complicated  slightly  by  the 
deformation  transformation  matrix  D.  The  product  of  the 
matrix  terms  involving  D  will  result  in  scaler  values, 
therefore  the  partial  derivatives  of  D  with  respect  to  each 
small  motion  generalized  coordinates  will  contribute 
separately.  This  requires  the  introduction  of  the  operator  5 
which  represents  a  column  vector  composed  of  the  matrix 
products,  where  the  derivatives  have  been  taken  with  respect 
to  each  small  motion  generalized  coordinate. 

The  symbol  DK  represents  the  deformation  transformation 
matrix  for  which  the  partial  derivative  of  matrix  D  have  been 
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taken  with  respect  to  each  small  motion  generalized 
coordinate  represented  by  K.  This  is  demonstrated  in 
equation  (C-15)  where  K  is  1.  Applying  the  operator  c  we  can 
express  the  scaler  terms  for  the  derivatives  of  D  for  each  K, 
1-6  by  equation  (C-16). 


0 

0                 0 

0 

d    =    2£    = 

Ul 
0 

0                 0 
0                 0 

0 
0 

0 

0                 0 

0 

p      T    T                i 

°  "eWDK=i 

T    T 
C   DXW^W   DK    = 

T    T 

D  we"DK=2 

■ 
■ 
■ 

T    T 

L  D  weWDK=6  J 

(C-15) 


(C-16) 


Separating  the   contributions  from   the  link  and  the  load 
we  present: 
For  the  link; 


9KE 


9u 


rf  =  J»[*TWTW*U  +   *TWTW(r+*U)]  dv 


glt^E]  =  Ju[*T(fiTw  +  wTfi)*5  +  ♦Vwfl"  ♦    *Vw*u 


aKE 
auT 


♦  ^(W1    W    +   &TW    )(r+*U)]    dv 


=    fy[*TWTW(r    +   *U)    +   *TWTW*U]    dv 


(C-17) 


(C-18) 


(C-19) 
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3-r|_— 7f?rJ    ~   — m    =   JUL*   w   w*u      +      *   W   W*U    +   *   W   W    (r+*U)j    dv 

(C-20) 
For    the    load: 

-TT  "  K&LDk","6cL*   ^DT8T"DkEl]      dvL  (C-21) 

*UK 

3UK 


r,D    (W      W    +   W1W)Dvr,+    r    D   W   WDvr, 1    dvT 
L  K- L       —  K~" LJ  L 


(C-22) 


3KE 

Tin  '  K&£DK8TflD|V  sW^tJ    dvL  (C-23' 

K 


hJ      Li         ^      L         f      f  T   T  Taj  t   t   T    ■■  T"T«T 

J±r ill    -  _  =    \u    Tr   D   W   WDr    +   r.D.1  ID   r    +   r   D   W   WD   r    + 

auK         auK 

E^TW      WDk^]    dvL  (C-24) 


Completing  our  development,  the  expressions  for  the 
potential  energy  are  presented  organized  as  the  previous 
material  and  followed  by  the  Lagrange  equations. 


PE    =  j        UTrTCrU   dx    +    fy(r    +   *U)TWTg_  dv    -   JuLr^DTWTg.  dv 


L 
(C-25) 


~|  =   "   JW  +   *U)TW*g_  dv    -   JwLr^DTwJg  dvL  (C-26) 
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Recalling  equations  (C-2) ,  (C-23),  and  (C-26)  we  can  make 
the  appropriate  substitutions  and  present  the  Lagrange 
equation  for  the  large  motion  generalized  coordinates. 

F   -  Jw[(r_+*U)T(WQW')(r+*U)  +  2(r_+*U)  Tw£w(*U)  + 

(r_+*U)Tw£w(*u')]  dv  +  JuL[(Dr^TWQW*(Dr^  + 
2(Dr£TWQ»(5r£  +  (Dr^W^CD'r^]  dvL  + 
-  Jw(r  +  *U)Tw£g.  dv  -  JuLr^DTwJg.  dvL  (C-27) 

We  can  simplify  this  expression  further  by  the 
substitution  of  the  integral  terms  defined  in  TABLE  II  and 
the  separation  of  the  second  time  derivative  of  the 
coordinate  transformation  matrix  W  by  its'  force 
contributions,  equation  (C-28).  The  contributions  to 
coriolis  and  centripetal  forces  are  represented  by  Wr,  termed 
the  residual  accelerations.  The  contribution  to  the  general 
forces  is  represented  by  We. 


W   =  W   +  W.  (C-28) 

r     w 

F     =  in(«X+  weV  +  Ii2(WX+  weVu  +  uTl2i(WeV  weV 
+uTl22(WJwr+  wJwe)U  +  2Ii2(Wefi)"  +  2UTi22(wJw)*u  +  iiaCi*mft' 

+    UTI22(WqW)u"  +  H5(DTwJwrD)  +  H5(DTWqW6D)  +  2H5(DTw£w6) 

+    H5(DTwJw5")  -  H^g.  -  UTH2w£g.  -  H^wjg.  (C-29) 
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Returning  to  the  development  of  the  small  motion  Lagrange 
equation  the   potential   energy   contribution  remains   to  be 
presented.    Again  as   in  the  previous  case  the  contributions 
of  the  point  load  and  the  link  will  be  developed  separately. 
For  the  link: 

2FJL.  =        rTCrU  dx   -    fy*TWTg_  dv  (C-30) 

dtSuTJ       auT      1? 


o  =     i22(wTw)5"  +  2i22(wTw)5  +  i21(«Twr)  +  i21(wTwQ)   + 

+  I22(WTWr)U  +  I22(WTWe)U  +  KXU  -  H2WTg_  (C-32) 


For  the  load: 
dPE 


auK 


=  -  JuLe£DkwT9-  dvL  (C-33) 


aTLTT-J       au-  +  au-  -  °  <c  34) 


auK        wwk 


T  Taj         t  T  ■ ■         T  T«  ■ 
0  =  Hc(Diw1WD)  +  HC(DAW1WD  )  +  HC(D  WAW  Dv)    + 

T    Taj  T    T 

H5(D1W1WDK)    -   H4DjJWa  (C-35) 

Collecting  terms  and  arranging   the  coefficients   the  two 
Lagrange  equations  can  be  written  as  the  equations  of  motion 
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for  the  large  and  small  generalized  coordinates.  It  is  these 
equations  which  must  be  solved  by  the  computer  simulation 
code. 

For  the  large  motion: 


Ii2(We"e)0  + 
uTl2i(WeV  + 
uTl22(WeVu  + 
H5  (DT»Xi» 


Ii2<"e")   + 

T  T 

0  I22(WeW)  + 

5   H5<DTWgWDK) 


0" 


F   "   IU<»eV    "   I12(WeWr,U    -   uTl21(M9Wr' 
2I12(Wgii)0    -   2UTI22(W^8)U    -   H5(DT»JwrD) 

H^a  +  uTH2w^a  +  H^w^a 


°Tl22(W9Wr)U 
2HC(D  W^WD)  + 


For  the  small  motion: 


(C-36) 


I22(wTwe)U  + 

c  H5(D^WTWQD) 


e 


I22(WAW)   + 
c  H5(D^WTWDK) 


2  I22(WTW)   + 
2  c  H5(D^WTWDK) 


■ 
U 


I22(W^r)  +  Kx  + 
c  H5(Dy«rDK) 


U 


[  -  I21(WTWr)  +  H2WTg.  +  c  H4D^WTg.  ] 


(C-37) 
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APPENDIX  D 
LISTING  OF  COMPUTER  SIMULATION  CODE 


******  ******  ***************** ******************************************* 

*  THIS  IS  THE  CONTROLLING  PROGRAM  FOR  BOOM.   IT  DEFINES  THE  ■ 
VARIABLES,  DETERMINES  THE  COEFFICIENTS  FOR  THE  EXPLICIT  * 

*  INTEGRATION,  AND  CALLS  THE  TWO  SEPARATE  INTEGRATOR  OF  THE  * 
LARGE  AND  SMALL  MOTION  EQUATIONS.  * 

*  * 

DEFINITION  OF  VARIABLES  * 

*  * 

A(2) H-D  JOINT  VARIABLE 

*  ALPHA(2) H-D  JOINT  VARIABLE  * 

AREA CROSS  SECTIONAL  AREA  OF  BOOM  * 

E YOUNG  '  S  MODULES  * 

EPSLON PARAMETER  FOR  INTEGRATION  COEFFICIENTS 

*  FINTIM FINISH  TIME  OF  INTEGRATION  * 

*  G SHEAR  MODULES  * 

GRAV GRAVITATIONAL  FORCE  (GLOBAL  COORDINATES)  * 

*  H TIME  STEP  USED  IN  THE  INTEGRATION  * 

IOTA PARAMETER  FOR  INTEGRATION  COEFFICIENTS 

IYY MOMENT  OF  INERTIA  ABOUT  THE  Y  AXIS 

122 MOMENT  OF  INERTIA  ABOUT  THE  2  AXIS  * 

MU(1) MASS  DENSITY  OF  BOOM  (KG/M**3)  * 

MU(2) POINT  MASS  AT  BOOM  TIP  (KG) 

Q LARGE  MOTION  GENERALI2ED  COORDINATE  9  * 

*  QDOT .LARGE  MOTION  ROTATIONAL  VELOCITY  * 

T INDEXING  VARIABLE  TO  PERMIT  TIMING  OF  PRINT   * 

*  TIME SEC  * 

*  Q2DOT LARGE  MOTION  ROTATIONAL  ACCELERATION  * 

*  TORK LARGE  MOTION  FORCING  FUNCTION  .  * 

*  U(l) SMALL  MOTION  AXIAL  DEFLECTION  * 

*  U(2) SMALL  MOTION  TRANSVERSE  DEFLECTION  (Y)  * 

*  U(3) SMALL  MOTION  TRANSVERSE  DEFLECTION  (2) 

U(4) SMALL  MOTION  TORSIONAL  DEFLECTION  * 

U(5) SMALL  MOTION  ROTATION  ABOUT  LOCAL  Y  AXIS 

U(6) SMALL  MOTION  ROTATION  ABOUT  LOCAL  2  AXIS      * 

*  * 

************************************************************************ 

INTEGER  PT,T 

REAL* 8  A ( 2 ) , D ( 2 ) , THETA ( 2 ) , ALPHA ( 2 ) , AREA, E, G , 
+        IYY, 122, MU(2), GRAV, 
+        U(6),UDOT(6),U2DOT(6), 
+      -  Q,QDOT,Q2DOT, 
+        EPSLON, IOTA, H, TIME, FINTIM, TORK 

COMMON  /BLOCK1/  A, D, THETA, ALPHA, AREA, E,G, IYY, 122, MU, GRAV 

COMMON  /BLK15/  Q, QDOT, Q2DOT 

COMMON  /BLK16/  U,UDOT,U2DOT 

GRAV=0 .DO 

TORK=100.D0 

*  MU(2)=44.1 
MU(2)=0.D0 
FINTIM  =2. DO 
H  =  .0006D0 
PT  =  34 

EPSLON  =  0. 25D0 
IOTA  =  0.5D0 
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A(1)=0.0 
D( 1)=0. 
THETA( 1)=Q 
ALPHA(1)=2.*ATAN(1.  ) 

A(2)=5.486 

D(2)=0. 

ALPHA(2)=0. 

AREA=1.9E-3 

E=70.E9 

G=26.0E9 

IYY=34.63E-6 

IZZ=34.63E-6 

MU(1)=2710.D0 

Q  =  O.DO 
QDOT  =  O.DO 
Q2DOT  =  O.DO 
DO  5  1=1,6 

U(I)=O.DO 

UDOT(I)=O.DO 
5       U2DOT(I)=O.DO 

THETA(2)=0.D0 
T=0 

DO  10  TIME=0,FINTIM,H 
T=T  +  1 

CALL  CCOEF(EPSLON,IOTA,H) 
CALL  INTSML(TIME,TORK) 
CALL  INTLRG(TIME.H,TORK) 
IF(TIME.EQ.0.DO)CALL  RINTl(TIME) 
IF(T.EQ.PT)THEN 

CALL  RINTHTIME) 
T=0 
ENDIF 
10   CONTINUE 

END 

SUBROUTINE  INTSMLCTIME, TORK) 
A*********************************************************************** 

»  * 

*  THIS  SUBROUTINE  IMPLICITLY  INTEGRATES  THE  SMALL  MOTION  EQUATION   * 

*  * 

*  ADDITIONAL  DEFINITION  OF  VARIABLES  * 

*  * 

*  Al-7 THE  COEFFICIENTS  FOR  NUMERICAL  INTEGRATION 

*  FN ...THE  SMALL  MOTION  FORCE  VECTOR 

*  FNN REFORMED  SMALL  MOTION  FORCE  VECTOR  * 

*  FQ THE  LARGE  MOTION  FORCE  VECTOR  * 

*  GN GYROSCOPIC  MATRIX  * 

*  KN STIFFNESS  MATRIX 

*  MN REFORMED  SMALL  MOTION  INERTIA  AND  COUPLING    * 

*  MATRIX  * 

*  MNN SMALL  MOTION  INERTIA  MATRIX 

*  MQN COUPLING  INERTA  MATRIX  * 

*  MQQ LARGE  MOTION  INERTIA  MATRIX  * 

*  TEMP1-5 TEMPORATY  MATRICES  FOR  INTERMEDIATE  VALUES  * 

*  PRD DUMMY  ARGUMENT  USED  IN  TROUBLE  SHOOTING  * 

*  WKAREA STORAGE  FOR  IMSL  LEQT2F  ROUTINE 

*  (USED  FOR  LINEAR  EQUATION  SOLVING)  * 

*  * 

a*********************************************************************** 


53 


REAL* 8  MQQ, MQN ( 6 ) , MNN ( 6 , 6 ) , GN ( 6 , 6 ) , KN ( 6 , 6 ) , FN ( 6 ) , FQ , 
+        U(6),UDOT(6),U2DOT(6), 
+        MN(6,6),FNN(6), 
+        A0,A1,A2,A3,A4,A5,A6,A7, 
+        WKAREA(54),PRD(6,6),TIME,TORK, 
+        TEMPI ( S ) , TEMP2 ( 6 ) , TEMP3 ( 6  )  , TEMP4 ( 6 ) , TEMP5 ( S . S ) 

COMMON  /BLK16/  U,UDOT,U2DOT 

COMMON  /BLOCK5/  MQQ, MQN, MNN, GN, KN,FN, FQ 

COMMON  /BLOCX7/  MN,FNN 

COMMON  /BLOCK8/  AO, Al, A2, A3, A4, A5, A6, A7 

CALL  TRANSF(TIME) 

CALL  WMATRX 

CALL  WDERIV 

IF( TIME. EQ. 0 . DO )CALL  LRGCOF( TORK ) 

CALL  SMLCOF 

CALL  REFORM  § 

DO  10  1=1,6 

TEMPKD     =   A0*U(D    +    A2*UDOT(I)     +    A3*U2DOT(I) 
10  TEMP2CI)     =    A1*U(I)     f    A4*UDOT(I)     +    A5*U2DOT(I) 

CALL  MATMUL(PRD,MN,6,6,TEMP1,6,1,TEMP3) 
CALL  MATMUL(PRD,GN,S,S,TEMP2,6,1,TEMP4) 

DO  20  1=1,6 

TEMPKD  =  FNN(I)  +  TEMP3U)  +  TEMP4  ( I ) 
TEMP2U)  =  U(I) 

DO  20  J=1.6 
20  TEMP5(I,J)  =  KN(I.J)  +  A0*MN(I,J)  +  A1*GN(I,J) 

M  =  1 
N  =  6 
IDGT  =  1 

CALL  LEQT2F ( TEMP5 , M, N. N , TEMPI , IDGT, WKAREA , IER ) 

DO    30    1=1,6 
30  U(I)    =    TEMPKD 

DO    40    1=1,6 

TEMP3(D     =    U2DOTCD 

U2DOT(D     =    AOMU(I)    -    TEMP2(D)    -    A2*UDOT(D    -    A3*TEMP3(I 
40  UDOT(I)     =    UDOT(I)     +    A6*TEMP3(D     +    A7*U2DOT(D 

END 
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SUBROUTINE  INTLRG  ( TIME, H, TORK) 
a*********************************************************************** 

*  * 

*  THIS  SUBROUTINE  INTEGRATES  THE  LARGE  MOTION  EQUATION  EXPLICITLY   * 

*  * 

*  ADDITIONAL  DEFINITION  OF  VARIABLES  * 

*  * 

*  Al-7 THE  COEFFICIENTS  FOR  NUMERICAL  INTEGRATION    * 

*  FN THE  SMALL  MOTION  FORCE  VECTOR  * 

*  FQ THE  LARGE  MOTION  FORCE  VECTOR  * 

*  FQN REFORMED  LARGE  MOTION  FORCE  * 

*  GN GYROSCOPIC  MATRIX  * 

*  KN STIFFNESS  MATRIX  * 

*  MN REFORMED  SMALL  MOTION  INERTIA  AND  COUPLING    * 

*  MATRIX  * 

*  MNN SMALL  MOTION  INERTIA  MATRIX  * 

*  MQN COUPLING  INERTA  MATRIX  * 

*  MQQ LARGE  MOTION  INERTIA  MATRIX  « 

*  TEMPI TEMPORATY  MATRIX  FOR  INTERMEDIATE  VALUES  * 

*  PRD ..DUMMY  ARGUMENT  USED  IN  TROUBLE  SHOOTING  * 

*  * 
*******»•»*«*******»***«-«*»»*«»*»*««**»**«*♦****«***«**.«***»«**,».,♦*. 

REAL* 8  MQQ, MQN ( 6 ) , MNN ( 6 , 6 ) , GN ( 6 , 6 ) , KN ( 6 , 6 ) , FN ( 6 ) , FQ, 

+  U(6),UDOT(6),U2DOT(6), 

+  Q,QDOT,Q2DOT, 

+  A0,A1,A2,A3,A4,A5,A6,A7, 

+  H, TEMPI, PRD(6,1), FQN, TORK, TIME 

COMMON  /BLOCK5/  MQQ, MQN, MNN, GN, KN, FN, FQ 

COMMON  /BLK15/  Q,QDOT,Q2DOT 

COMMON  /BLK16/  U,UDOT,U2DOT 

COMMON  /BLOCK8/  AO , Al, A2, A3 , A4, A5 , A6, A7 

Q  =  Q  +  H*QDOT  +  A3*Q2DOT/A0 
QDOT  =  QDOT  +  A6*Q2DOT 

CALL  TRANSF(TIME) 
CALL  WMATRX 
CALL  WDERIV 
CALL  LRGCOF(TORK) 

v 

TEMP1=0.D0 
DO  10  1=1,6 
10      TEMPI  =  TEMPI  +  MQN(I)*U2DOT(I) 

FQN  =  FQ  -  TEMPI 

Q2DOT  =  ( l.DO/MQQ)*FQN 
QDOT  =  QDOT  +  Q2DOT*A7 
Q  =  Q  +  Q2DOT/A0 
END 


55 


SUBROUTINE  LRGCOF(TORK) 
a*********************************************************************** 

*  * 

*  THIS  SUBROUTINE  CALCULATES  VALUES  FOR   COEFFICIENT  MATRICES  FOR  * 

*  LINK  ONE.  * 

*  * 

*  ADDITIONAL  DEFINITION  OF  VARIABLES  * 

*  * 

*  Al-7 THE  COEFFICIENTS  FOR  NUMERICAL  INTEGRATION  * 

*  CONST SINGLE  VALUED  TERM  USED  IN  SUMMATION  FOR 

*  FQ  VALUE  * 

*  DC DERIVATIVE  OF  DFORM  MATRIX  WRT  SMALL  MOTION  * 

*  COORDINATE  OF  INTEREST  * 

*  DCT TRANSPOSE  OF  DC  MATRIX  * 

*  DDOT THE  RATE  OF  DEFORMATION  TRANSFORMATION  MATRIX* 

*  DDOTT THE  TRANSPOSE  OF  THE  DDOT  MATRIX  * 

*  DFORM THE  DEFORMATION  TRANSFORMATION  MATRIX  * 

*  DFORMT THE  TRANSPOSE  OF  THE  DFORM  * 

*  FN THE  SMALL  MOTION  FORCE  VECTOR  * 

*  ,     FQ THE  LARGE  MOTION  FORCE  VECTOR  * 

*  FQN REFORMED  LARGE  MOTION  FORCE  * 

*  GN GYROSCOPIC  MATRIX  * 

*  KN STIFFNESS  MATRIX  * 

*  MAT*** VARIOUS  DIMENSIONED  MATRICES  OF  INTERMEDIATE* 

*  VALUES  (***)  FOR  THE  VARIOUS  SIZES  * 

*  MN REFORMED  SMALL  MOTION  INERTIA  AND  COUPLING  * 

*  MATRIX  * 

*  MNN SMALL  MOTION  INERTIA  MATRIX  * 

*  MQN COUPLING  INERTA  MATRIX  * 

*  MQQ LARGE  MOTION  INERTIA  MATRIX  * 

*  TEMPI : TEMPORATY  MATRIX  FOR  INTERMEDIATE  VALUES  '* 

*  ONE ARRAY  OF  UNITY  NECESSARY  IN  OBTAINING 

*  DERIVATIVE  OF  DFRM  MATRIX  .* 

*  PRD DUMMY  ARGUMENT  USED  IN  TROUBLE  SHOOTING  * 

*  W HARTENBERG  DENAVIT  TRANSFORMATION  MATRIX  * 

*  WPRM DERIVATIVE  OF  W  WRT  9  * 

*  WDOT DERIVATIVE  OF  W  WRT  TIME  * 

*  WRESID RESIDUAL  ACCELERATION  TRANSFORMATION  * 

*  Wl-14 VARIOUS  KERNELS  FOR  THE  BOOM  VOLUME  INTEGRAL  * 

*  XL LOWER  LIMIT  OF  BOOM  VOLUME  INTEGRATION  * 

*  XU UPPER  LIMIT  OF  BOOM  VOLUME  INTEGRATION  * 

*  Y,Z DUMMY  ARGUMENTS  * 

*  * 

INTEGER  INDEX, NR, NC, NR2, NC2, TERM 

REAL* 8  A  (  2  ) ,  D  (  2  ) , THETA ( 2 ) , ALPHA (  2 ) , AREA , E , G , 

+  IYY,IZZ,MU(2),GRAV, 

+  U(6),UDOT(6),U2DOT(6), 

+  UT(6), 

+  Q,QDOT.Q2DOT. 

+  W(4,4), 

+  WPRM(4,4),WDOT(4,4),WRESID(4,4), 

+  WD1(4,4),WD2(4,4),WD3(4,4),WD4(4,4),WD6(4,4), 

+  WD9(4,4).WD10(4,4).WD13(4,4),WD14(4,4), 

+  MQQ , MQN (  6 )  ,  MNN ( 6 , 6 ) , GN ( 6 , 6 ) , KN ( 6 , 6  ) ,  FN  (  6  )  ,  FQ , 

+  MAT(6),MAT4A(4),MAT4B(4),MATS4A(6.4),MAT64B(6,4) , MATPRD( 6, 6 ) , 

+  CONST,Z(6,6),Y(6,6),XL,XU,ONE(6),MAT44A(4,4),MAT44B(4,4), 

+  PTMASS(4),GRAVEC(4),WT(4,4),WPRMT(4,4),PRD(6,6) , DC (4,4), 

+  DFRM(4,4),DFRMT(4,4).DDOT(4,4),DDOTT(4,4),DCT(4,4),TORK 


56 


COMMON  /BLOCK1/  A,D, THETA.  ALPHA, AREA, E,G, I YY, IZZ,MU,GRAV 

COMMON  /BLK15/  Q,QDOT,Q2DOT 

COMMON  /BLK16/  U,UDOT,U2DOT 

COMMON  /BLOCK 3/  W 

COMMON  /BLOCK4/  WD1, WD2. WD3, WD4, WD6, WD9, MD10, WD13, WD14 

COMMON  /BLOCK5/  MQQ,MQN, MNN,GN, KN, FN,FQ 

COMMON  /BLOCK6/  WPRM, WDOT, WRESID 

XL=-A(2) 
XU=0.D0 

MQQ=0.D0 
FQ=TORK 
DO  3  1=1,6 
3       MQN(I)=0.D0 

DO  5  J =1,4 

GRAVEC(J)=0.D0 
5       PTMASS(J)=0.D0 

GRAVEC(4)=GRAV 

PTMASS(1)=MU(2) 

TERM=0 

DO  7  J=l,6 
7       ONE(J)=l.D0 

CALL  DFORM(l,7,U,DFRM) 

CALL  DFORM(0,7,UDOT,DDOT) 

CALL  TRANSP(TERM,DFRM,4,4,DFRMT) 

CALL  TRANSP(TERM,DDOT,4,4,DDOTT) 

'j 

-,  CALL  TRANSP(TERM,U,6, 1,UT) 
CALL  TRANSP(TERM,W,4,4,WT) 
CALL  TRANSP(TERM,WPRM,4,4,WPRMT) 

a******************************************************* ***••*••***•*** 

*  CREATING  THE  COEFFICIENT  MATRIX  MQQ 
a**********-************************************************************ 

INDEX  =  11 

CALL  BINT( TERM, XL, XU,Z,Y. 1,1, INDEX. WD6, CONST) 

MQQ=MQQ+CONST 

INDEX  =  12 

CALL  BINT ( TERM, XL, XU,Z,Y, 1,6, INDEX, WD6, MAT) 

CALL  MATMUH  PRD, MAT, 1 , 6 , U , 6 , 1 , CONST ) 

MQQ=MQQ+CONST 

INDEX  =  21 

CALL  BINT( TERM, XL, XU,Z,Y, 6,1, INDEX, WD6, MAT) 

CALL  MATMUH PRD, UT, 1,6, MAT, 6,1, CONST) 

MQQ=MQQ+CONST 

INDEX  =  22 

CALL  BINT ( TERM, XL, XU, Z, Y, 6, 6, INDEX, WD6 , MATPRD ) 

CALL  MATMUL(PRD,MATPRD,6,6,U,6, 1,MAT) 

CALL  MATMUL ( PRD , UT, 1,6, MAT, 6,1, CONST ) 

MQQ=MQQ+CONST 

*  THE  CONTRIBUTION  OF  THE  POINT  LOAD 

CALL  MATMUL(PRD,WD6,4,4,DFRM,4,4,MAT44A) 
CALL  MATMUL(PRD,DFRMT,4,4,MAT44A,4,4,MAT44B) 

*  ONLY  ELEMENT  MAT44B(1,1)  MILL  SURVIVE  FOLLOWING  POST  MULTIPLICATION  BY 

*  LOCAL  POSITION  VECTOR  AND  PRE  MULTIPLICATION  BY  ITS  TRANSPOSE. 
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CONST=MU(2)*MAT4  4B(l,l) 
MQQ=MQQ+CONST 

•a********************************************************************* 

*  CREATING  THE  COEFFICIENT  MATRIX  HQN 
A********************************************************************** 

INDEX=12 

CALL  BINT ( TERM, XL, XU,Z,Y. 1,6, INDEX, WD2, MAT) 
DO  8  1=1,6 

8  MQN(I)=MQN(I)  +  MAT(I) 

INDEX  =  22 

CALL  BINT (TERM, XL, XU,Z,Y, 6, 6, INDEX, WD2, MATPRD) 

CALL  MATMUL(PRD,UT,1,6,MATPRD,6,6,MAT)- 

DO  9  1=1,6 

9  MQN(I)  =  MQN(I)  +  MAT(I) 

*  THE  CONTRIBUTION  OF  THE  POINT  LOAD 

DO  10  K=l,6 

CALL  DFORM(0,K,ONE,DC) 

CALL  MATMUL(PRD,WD2,4,4.DC,4,4,MAT44A) 

CALL  MATMUL(PRD,DFRMT,4,4,MAT44A,4,4,MAT44B) 

*  ONLY  ELEMENT  MAT44B(1,1)  MILL  SURVIVE  FOLLOWING  POST  MULTIPLICATION  BY 

*  LOCAL  POSITION  VECTOR  AND  PRE  MULTIPLICATION  BY  ITS  TRANSPOSE. 

10  MQN(K)=MQN(K)  +  MU ( 2 ) *MAT44B ( 1, 1 ) 


*  CREATING  THE  COEFFICIENT  MATRIX  FQ 

INDEX  =  11 

CALL  BINT(TERM, XL, XU,Z,Y, 1,1, INDEX, WD14, CONST) 

FQ=FQ-CONST 

INDEX  =  12 

CALL  BINT ( TERM, XL, XU,Z,Y, 1,6, INDEX. WD10, MAT) 
CALL  MATMUL(PRD,MAT,l,6,UDOT,6, 1, CONST) 
FQ=FQ-2. *CONST 

CALL  BINT ( TERM, XL, XU.Z.Y, 1,6, INDEX, HD14, MAT) 
CALL  MATMUL(PRD,MAT,l,6,U,6,l,CONST) 
FQ=FQ-CONST 

INDEX  =  21 

CALL  BINT ( TERM, XL,  XU,Z,Y, 6,1, INDEX, WD1 4,  MAT) 

CALL  MATMUL(PRD,UT, 1,6, MAT, 6,1, CONST) 

FQ=FQ-CONST 

INDEX  =  22 

CALL  BINT(TERM,XL,XU,Z,Y,6,6,INDEX,WD10,MATPRD) 

CALL  MATMUL(PRD,MATPRD,6,6,UDOT,6,l,MAT) 

CALL  MATMUL(PRD,UT, 1,6, MAT, 6,1, CONST) 

FQ=FQ-2. *CONST 

CALL  BINT ( TERM, XL, XU, Z, Y, 6, 6, INDEX, WD14 , MATPRD ) 
CALL  MATMULCPRD, MATPRD, 6, 6,U, 6,1, MAT) 
CALL  MATMULCPRD, UT, 1,6, MAT, 6,1, CONST) 
FQ=FQ-CONST 

*  THE  CONTRIBUTION  OF  THE  POINT  LOAD 
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CALL  HATMUL(PRD,WD14,4,4,DFRM,4,4,MAT44A) 
CALL  MATMUL(PRD,DFRMT,4,4,MAT44A,4,4,MAT44B) 

*  ONLY  ELEMENT  MAT44B(1,1)  WILL  SURVIVE  FOLLOWING  POST  MULTIPLICATION  BY 

*  LOCAL  POSITION  VECTOR  AND  PRE  MULTIPLICATION  BY  ITS  TRANSPOSE. 

CONST=MU(2)*MAT44B(l,l) 
FQ=FQ-CONST 

CALL  MATMUL(PRD,WD10,4,4,DDOT,4,4,MAT44A) 
CALL  MATMUL(PRD,DFRMT,4,4,MAT44A,4,4,MAT44B) 

*  ONLY  ELEMENT  MAT44B(1,1)  WILL  SURVIVE  FOLLOWING  POST  MULTIPLICATION  BY 

*  LOCAL  POSITION  VECTOR  AND  PRE  MULTIPLICATION  BY  ITS  TRANSPOSE. 

CONST=2 . DO  *MU ( 2 ) *MAT44B  < 1,1) 
FQ=FQ-CONST 

*  THE  CONTRIBUTION  OF  POTENTIAL  ENERGY 

INDEX  =  10 

*  HI,  WD1  WILL  BE  CARRIED  AS  A  DUMMY  ARGUMENT 

CALL  BINT ( TERM, XL, XU,Z,Y, 1,4, INDEX, WD1,MAT4 A) 
CALL  MATMUL(PRD,MAT4A, 1, 4, WPRMT, 4, 4, MAT4B ) 
CALL  MATMUL(PRD,MAT4B,1,4,GRAVEC,4,  1, CONST) 
FQ=FQ+CONST 

INDEX  =  20 

*  H2,  WD1  WILL  BE  CARRIED  AS  A  DUMMY  ARGUMENT 

CALL  BINT(TERM,XL,XU,Z,Y,6,4,INDEX,WD1,MAT64A) 
CALL  MATMUL(PRD,UT,1,6,MAT64A,6,4,MAT4A) 
CALL  MATMUL(PRD,MAT4A, 1, 4, WPRMT, 4, 4,MAT4B) 
CALL  MATMUL(PRD,MAT4B,l,4,GRAVEC,4,l,CONST) 
FQ=FQ+CONST 

*  THE  CONTRIBUTION  OF  THE  POINT  LOAD  TO  THE  POTENTIAL  ENERGY 

CALL  MATMUL(PRD,DFRMT,4,4,WPRMT,4,4,MAT44A) 
CALL  MATMUH  PRD, MAT44A, 4, 4, GRAVEC, 4/ 1, MAT4B) 
CALL  MATMUL(PRD,PTMASS,l,4,MAT4B,4,l,CONST) 
FQ=FQ+CONST 

RETURN 
END 
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SUBROUTINE  SHLCOF 

A*********************************************************************** 

*  » 

*  THIS  SUBROUTINE  CALCULATES  VALUES  FOR   COEFFICIENT  MATRICES  FOR     * 

*  LINK  ONE.  * 

*  * 

*  ADDITIONAL  DEFINITION  OF  VARIABLES  * 

*  * 

*  Al-7 THE  COEFFICIENTS  FOR  NUMERICAL  INTEGRATION    * 

*  CONST SINGLE  VALUED  TERM  USED  IN  SUMMATION  FOR      * 

*  FQ  VALUE  * 

*  DC DERIVATIVE  OF  DFORM  HATRIX  WRT  SMALL  MOTION   * 

*  COORDINATE  OF  INTEREST  * 

*  DCT TRANSPOSE  OF  DC  MATRIX  * 

*  DDOT THE  RATE  OF  DEFORMATION  TRANSFORMATION  MATRIX* 

*  DDOTT THE  TRANSPOSE  OF  THE  DDOT  MATRIX  * 

*  DFORM THE  DEFORMATION  TRANSFORMATION  MATRIX  * 

*  DFORMT THE  TRANSPOSE  OF  THE  DFORM  * 

*  FN THE  SMALL  MOTION  FORCE  VECTOR  * 

*  FQ THE  LARGE  MOTION  FORCE  VECTOR  * 

*  FQN REFORMED  LARGE  MOTION  FORCE  * 

*  GN GYROSCOPIC  MATRIX  * 

*  KN STIFFNESS  MATRIX  * 

*  MAT*** VARIOUS  DIMENSIONED  MATRICES  OF  INTERMEDIATE  * 

*  VALUES  (***)  FOR  THE  VARIOUS  SIZES  * 

*  MN REFORMED  SMALL  MOTION  INERTIA  AND  COUPLING    * 

*  MATRIX  * 

*  MNN SMALL  MOTION  INERTIA  MATRIX  * 

*  MQN COUPLING  INERTA  MATRIX  * 

*  MQQ LARGE  MOTION  INERTIA  MATRIX  * 

*  TEHP1 TEMPORATY  MATRIX  FOR  INTERMEDIATE  VALUES  * 

*  ONE ARRAY  OF  UNITY  NECESSARY  IN  OBTAINING  * 

*  DERIVATIVE  OF  DFRM  MATRIX  * 

*  PRD DUMMY  ARGUMENT  USED  IN  TROUBLE  SHOOTING  * 

*  W HARTENBERG  DENAVIT  TRANSFORMATION  MATRIX  * 

*  WPRM -.DERIVATIVE  OF  W  WRT  9  * 

*  WDOT DERIVATIVE  OF  W  WRT  TIME  * 

*  WRESID RESIDUAL  ACCELERATION  TRANSFORMATION  * 

*  Wl-14 VARIOUS  KERNELS  FOR  THE  BOOM  VOLUME  INTEGRAL  * 

*  XL LOWER  LIMIT  OF  BOOM  VOLUME  INTEGRATION  * 

*  XU UPPER  LIMIT  OF  BOOM  VOLUME  INTEGRATION  * 

*  Y,Z DUMMY  ARGUMENTS  * 

*  * 

INTEGER  INDEX, NR,NC,NR2^ NC2, TERM 

REAL*8  A ( 2 ) , D( 2 ) , THETA ( 2 ) , ALPHA ( 2 ) , AREA, E,G. 
+        IYY,IZZ,MU(2),GRAV, 
+        W(4,4), 

+        U(6),UDOT(6),U2DOT(6), 
+        UT(6), 
+        Q,QDOT,Q2DOT, 

+        WPRM(4,4),WDOT(4,4),WRESID(4,4), 
+        WD1(4,4),WD2(4,4),WD3(4,4),WD4(4,4),WD6(4,4), 
+        WD9(4,4),WD10(4,4),WD13(4,4),WD14(4,4), 
*        MQQ,MQN(6),HNN(6,6),GN(6,6),KN(6,6) ,FN(6),FQ, 

+  MAT(6),MAT4A(4).HAT4B(4),MAT64A(6,4),MAT64B(6,4).MATPRD(6,6), 
+        CONST,Z(6,6),Y(6,6),XL,XU,ONE(6),MAT44A(4,4),MAT44B(4,4), 
+        PTHASS(4),GRAVEC(4),WT(4,4),WPRHT(4,4).PRD(6,6),DC(4,4), 
+        DFRH(4,4),DFRMT(4,4),DOOT(4,4),DDOTT(4,4),DCT(4,4) 

COMMON  /BLOCK1/  A, D, THETA, ALPHA, AREA, E,G, I YY, IZZ, MU,GRAV 

COMHON  /BLK15/  Q,QDOT,Q2DOT 

COHMON  /BLK16/  U,UDOT,U2DOT 

COMHON  /BLOCK3/  W 

COMMON  /BLOCK4/  WD1, WD2, WD3, WD4, WD6, WD9, WD10, WD13, WD14 

COMMON  /BLOCK5/  MQQ, MQN, MNN,GN, KN, FN, FQ 

COMMON  /BLOCK6/  WPRM, WDOT, WRESID 
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XL=-A(2) 
XU-O.DO 
DO  3  1-1,6 
FN(I)=O.DO 
DO  3  J=l,6 

GN(I.J)=O.DO 
KN(I,J)=»O.DO 
3  MNN(I,J)=0.D0 

DO  5  J»l,4 

GRAVEC(J)=O.DO 
5       PTHASS(J)=0.D0 
GRAVEC(4)=GRAV 
PTMASS(1)=MU(2) 
TERM=0 
DO  7  J=l,6 
7       ONE(J)=l.D0 

CALL  DFORM(l,7,U,DFRM) 

CALL  DFORM(0,7,UDOT,DDOT) 

CALL  TRANSP(TERM,DFRH,4,4,DFRMT) 

CALL  TRANSP(TERM,DDOT,4,4,DDOTT)   • 

CALL  TRANSP(TERM.U,6,1,UT) 

CALL  TRANSP(TERM,W,4,4,WT) 

CALL  TRANSP(TERN,WPRM,4,4,WPRMT) 

*  CREATING  THE  COEFFICIENT  MATRIX  MNN 
*»***»***»«*»»«•««»««**»««*««««««««*»»»»»»»*»»»,», ********************** 

INDEX  =«  22 

CALL  BINT(TERH,XL,XU,Z,Y,6,6,INDEX,WD1,MATPRD) 
DO  12  1-1,6 

DO  12  J=l,6 
12  HNN(I,J)-MNN(I,J)  +  MATPRD(I,J) 

*  THE  CONTRIBUTION  OF  THE  POINT  LOAD 

DO  15  K=-l,6 

CALL  DFORM(0,K,ONE,DC) 

CALL  TRANSP(TERH,DC,4,4,DCT) 

CALL  MATHUL(PRD,DCT,4,4,WD1,4,4,MAT44A) 

CALL  DFORH(0,K.ONE,DC) 

CALL  MATMUL(PRD,MAT44A,4,4,DC,4,4,MAT44B) 

*  ONLY  ELEMENT  HAT44B(1,1)  MILL  SURVIVE  FOLLOWING  POST  MULTIPLICATION  BY 

*  LOCAL  POSITION  VECTOR  AND  PRE  MULTIPLICATION  BY  ITS  TRANSPOSE. 

15      HNN(K,K)=MNN(K,K)  +  HU( 2) »MAT44B< 1, 1) 

•a********************************************************************* 

*  CREATING  THE  COEFFICIENT  MATRIX  GN      (UDOT) 

a********************************************************************** 

INDEX  =  22 

CALL  BINT(TERH,XL,XU,Z,Y,6,6,INDEX,WD9,MATPRD) 

DO  22  1-1,6 

DO  22  J=l,6 
22  GN(I,J)=GN(I,J)+2*MATPRD(I,J) 

*  THE  CONTRIBUTION  FROM  THE  POINT  LOAD 

DO  30  J»l,6 

CALL  DFORM(0,J,ONE.DC) 

CALL  MATMUL(PRD,WD9,4,4,DC,4,4,MAT44A) 

CALL  TRANSP(TERM,DC,4,4,DCT) 

CALL  MATMUL(PRD,DCT,4,4,MAT44A,4,4,MAT44B) 

*  ONLY  ELEMENT  MAT44B(1,1)  WILL  SURVIVE  FOLLOWING  POST  MULTIPLICATION  BY 

*  LOCAL  POSITION  VECTOR  AND  PRE  MULTIPLICATION  BY  ITS  TRANSPOSE. 
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C0NST=2*MU( 2 ) *MAT44B( 1, 1) 
30      GN(J,J)=CN(J,J)+CONST 


••••••••a********************************************************** 


a********************************************************************** 

*  CREATING  THE  COEFFICIENT  MATRIX  KN(U) 

»«  »» 

INDEX  =»  22 

CALL  BINT(TERH,XL,XU,Z,Y,6,6,INDEX.HD13,MATPRD) 

DO  24  1=1,6 
DO  24  J=l,6 
24  KN(I,J)=»KN(I,J)+MATPRD(I,J) 

*  THE  CONTRIBUTION  OF  POTENTIAL  ENERGY  FROH  THE  LINK 

*  Kl,  HDl  HILL  BE  CARRIED  AS  A  DUMHY  ARGUMENT 

INDEX  =  1 

CALL  BINT(TERM,XL,XU,Z,Y,6,6,INDEX,HD1,MATPRD) 

DO  28  1-1,6 

DO  28  J  =  l,6 
28  KN(I,J)=KN(I,J)+HATPRD(I,J) 

*  THE  CONTRIBUTION  FROH  THE°  POINT  LOAD 

DO  32  J=l,6 

CALL  DFORM(0,J,ONE.DC)  • 

CALL  MATMUL(PRD,HD13, 4, 4,DC, 4, 4,MAT44A) 

CALL  DFORH(0,J,ONE,DC) 

CALL  TRANSP(TERM,OC,4,4,DCT) 

CALL  HATHUL(PRD,DCT,4.4,MAT44A,4,4,HAT44B) 

*  ONLY  ELEMENT  MAT44B(1,1)  HILL  SURVIVE  FOLLOHING  POST  MULTIPLICATION  BY 

*  LOCAL  POSITION  VECTOR  AND  PRE  MULTIPLICATION  BY  ITS  TRANSPOSE. 

CONST=MU( 2 ) »HAT44B( 1, 1 ) 
32      KN(J,J)=KN(J,J)+CONST 

*  CREATING  THE  COEFFICIENT  MATRIX  FN 

INDEX  =  21 

CALL  BINT (TERM, XL, XU„Z.Y,S,1, INDEX, HD13, MAT) 
DO  20  J=l,6 
20      FN(J)=FN(J)-MAT(J) 

*  THE  CONTRIBUTION  OF  POTENTIAL  ENERGY  FROH  THE  LINK 

*  H2,  HD1  HILL  BE  CARRIED  AS  A  DUHMY  ARGUMENT 

INDEX  =  20 

CALL  BINT ( TERM, XL, XU, Z, Y, 6, 4 , INDEX, HD1 . MATS4A ) 
CALL  MATMUL(PRD,MAT64A,6,4,HT,4,4,HAT64B) 
CALL  MATMUL(PRD,MAT64B,6.4,GRAVEC,4, 1,MAT) 
DO  26  J=l,6 
26      FN(J)=FN(J)+MAT(J) 

*  THE  CONTRIBUTION  OF  THE  POINT  LOAD  TO  THE  POTENTIAL  ENERGY 

DO  36  J=l,6 

CALL  DFORM(0.J,ONE,DC) 
CALL  TRANSP(TERM,DC,4,4,DCT) 
CALL  MATMUL(PRD,DCT,4,4,HT,4,4,MAT44A) 
CALL  MATHUL(PRD,MAT44A,4»4,GRAVEC,4,1,MAT4A) 
CALL  HATMUL(PRD,PTHASS, 1,4,MAT4A,4, 1, CONST) 
36      FN(J)=FN(J)+CONST 

RETURN 
END 
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SUBROUTINE  WMATRX 

ft*********************************************************************** 

*  * 

*  THIS  SUBROUTINE  MULTIPLIES  THE  HARTENBERG-DENAVIT  HATRICES        * 

*  TOGETHER  TO  OBTAIN  THE  COORDINATE  TRANSFORMATION  MATRIX  * 

*  • 

*  ADDITIONAL  DEFINITION  OF  VARIABLES  * 

*  * 

*  HD THE  HARTENBERG  DENAVIT  TRANSFORMATION  MATRIX  * 

*  W TRANSFORMATION  MATRIX  * 

*  * 

a*********************************************************************** 

REAL*8  W(4,4),HD(2,4,4),WSUM(4,4),SUM 

COMMON  /BLOCX2/  HD 
COMMON  /BLOCK3/  W 

DO  10  L=l,4 

DO  10  M  =  l,4 

WSUM(L,H)=0.D0 
10  W(L,M)=0.D0 

DO  40  J=l,4 

DO  30  K  =  l,4 
SUM=0.0 
DO  20  H  =  l,4 

SUM  =  SUM  +  HD(1,J,M)  *  HD(2,M,K) 
20  CONTINUE 

WSUH(J,K)=SUH 
30      CONTINUE 
40   CONTINUE 

DO  50  1-1,4 
DO  50  J=l,4 
50  W(I,J)=WSUH(I,J) 

END 

SUBROUTINE  WDERIV 

*  * 

*  THIS  SUBROUTINE  CREATES  THE  DERIVATIVES  OF  THE  COORDINATE  MATRIX  * 
»  W  WRT  Q  THE  LARGE  MOTION  GENERALIZED  COORDINATE  9  WPRM  * 

*  AND  TIME  WDOT  AND  WRESID.  * 

»  * 

*  ADDITIONAL  DEFINITION  OF  VARIABLES  * 

»  * 

*  HD THE  HARTENBERG  DENAVIT  TRANSFORMATION  MATRIX  * 

*  HDPRM DERIVATIVE  OF  HD(l)  WRT  9  * 

*  W TRANSFORMATION  MATRIX  • 

*  WT TRANSPOSE  OF  THE  TRANSFORMATION  HATRIX     * 

*  WPRM DERIVATIVE  OF  W  WRT  9  * 

*  WPRMT TRANSPOSE  OF  THE  DERIVATIVE  OF  W  WRT  9  * 

*  WDOT DERIVATIVE  OF  W  WRT  TIME  * 

*  WDOTT TRANSPOSE  OF  THE  DERIVATIVE  OF  W  WRT  TIME  * 

*  WRESID RESIDUAL  ACCELERATION  TRANSFORMATION  * 

*  WRSIDT TRANSPOSE  OF  THE  RESIDUAL  ACCELERATION  * 

*  TRANSFORMATION  * 

*  Wl-14 VARIOUS  KERNELS  FOR  THE  BOOM  VOLUME  INTEGRAL  * 

*  * 
a*********************************************************************** 

INTEGER  TERM 

REAL*8  W(4,4),HD(2,4,4),WPRM(4,4),WDOT(4,4),WRESID(4,4), 
+        Q,QDOT,Q2DOT, 

+        WT( 4, 4 ) , WPRMT ( 4, 4 ) , WDOTT ( 4, 4 ) , WRSIDT( 4, 4 ) , 
+        HDPRM(4,4),HD2(4,4),PRD<4,4), 

+        WD1(4,4),WD2(4.4),WD3(4,4),WD4(4,4),WD6(4,4),WD9(4,4), 
+        WD10(4,4),WD13(4,4),WD14(4,4) 
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COMMON  /BLK15/  Q,QDOT,Q2DOT 

COMMON  /BLOCK2/  HD 

COMMON  /BLOCK3/  U 

COMMON  /BLOCK4/  WD1, WD2, WD3, WD4 , WD6, WD9, WD10, WD13, WD14 

COMMON  /BLOCK6/  WPRM, WDOT, WRESID 

DO  5  K=l,4 
DO  5  L=l,4 

HDPRM(K,L)=0.0 
HD2(K,L)=HD(2,K,L) 
WDOT(K,L)=0.0 
5  WRESID(K,L)=0.0 

DO  10  K  =  1,4 

HDPRM(2,K)=-HD(1,3,K) 
10  HDPRH(3,K)=  HD(1,2,K) 

CALL  MATHUL(PRD,HDPRM,4,4,HD2,4,4,WPRM) 

*   THE  FOLLOWING  MATRIX  HAS  A  QDOT  COEFFICIENT 

DO  15  1-1,4 

DO  15  J=l,4 

IF(ABS(WPRM( I, J)). LE.l.E-30)  THEN 

WDOT(I,J)=0.D0 
ELSE 

WDOT ( I , J ) =WPRH ( I , J ) *QDOT 
ENDIF 
15      CONTINUE 

*  CREATING  THE  RESIDUAL  ACCELERATION  TRANSFORMATION  MATRIX 

*  THIS  HATRIX  HAS  A  QDOT**2  COEFFICIENT 

DO  20  K  >  1,4 

IF(ABS(WDOT(2,K) ) .LE.l.E-30)  THEN 

WRESID(2,K)=O.D0 
ELSE 

WRESID(2,K)=-WDOT(3,K)*QDOT 
ENDIF 
IF(ABS(WDOT(3,K) ) .LE.l.E-30)  THEN 

WRESID(3,K)=0.D0 
ELSE 

WRESID(3,K)=WDOT(2,K)*QDOT 
ENDIF 
20      CONTINUE 

CALL  TRANSP(TERM,W,4,4,WT) 
CALL  TRANSP(TERM,WPRM,4,4,WPRMT) 
CALL  TRANSP(TERM,WDOT,4,4,WDOTT) 
CALL  TRANSP(TERH,WRESID,4,4,WRSIDT) 

CALL  MATMUL( PRD, WT, 4.4,W,4,4,WD1) 
CALL  MATMULC  PRD . WPRMT, 4 , 4 , W , 4 , 4 , WD2 ) 
CALL  MATMUL(PRD,WDOTT,4,4,W,4,4,WD3) 
CALL  HATMUL(PRD,WRSIDT,4,4.W,4,4,WD4) 

CALL  MATMUL(PRD,WPRHT,4,4,WPRM,4,4,WD6) 

CALL  MATMUL(PRD,WT,4,4,WDOT,4,4,WD9) 
CALL  HATMUL(PRD,WPRMT,4,4,WDOT,4,4,WD10) 

CALL  MATMUL(PRD,WT,4,4,WRESID,4,4,WD13) 
CALL  MATMUL(PRD, WPRMT, 4, 4, WRESID, 4, 4, WD14 ) 

END 
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SUBROUTINE  INTPRD  (TERM, INDEX, WD, X, ROW, COL, PROD) 

*  * 

*  THIS  SUBROUTINE  CREATES  THE  PRODUCT  ARRAY  OF  THE  SHAPE  FUNCTIONS  * 

*  AND  THE  COORDINATE  TRANSFORMATION  MATRICES  W(I)  * 

* 
•a********************************************************************* 

INTEGER  INDEX, ROW, COL, TERM 

REAL*8  PROD(ROW,COL),R(4),P(4,6),WD(4,4), 
+        A ( 2 ) , D ( 2 ) , THETA ( 2 ) , ALPHA ( 2 ) , AREA, E, G, 
+        IYY,IZZ,MU(2),GRAV, 

+     P1(4,6),P1T(6,4),GMAT(4,6),GT(6,4),G1(4,6),G1T(6,4), 
+      WMAT(4,4),C(4,4),RC(4,1), 
+      PRD(6,6),X 

COMMON  /BLOCX1/  A, D, THETA, ALPHA, AREA, E,G, IYY, IZZ,MU,GRAV 

IF  (INDEX  .EQ.  11)  THEN 

CALL  VECTR(X  R) 

PROD(l,l)=WD(l,l)  +  R(2)*WD(2,1)  +  R(2)*WD(1,2) 

PROD(l,l)=PROD(l,l)+R(2)«*2»WD(2,2)+IYY*WD(3,3)/(MU(l)*AREA) 

PROD(l,l)«PROD(l,l)+IZZ*WD(4,4)/(HU(l)*AREA) 
END  IF 

IF  (INDEX  .EQ.  12)  THEN 
CALL  VECTR(X,R) 
CALL  C0EFP1(X,P1) 

PR0D(1,1)=WD(1,2)»P1(2,1)+R(2)*WD(2,2)*P1(2,1) 
PR0D(1,2)=WD(1,3)*P1(3,2)+R(2)*WD(2,3)*P1(3,2) 
PR0D(1,3)=WD(1,4)*P1(4,3)+R(2)*WD(2,4)*P1(4,3) 
PR0D(1,4)=(-IZZ*WD(4,3)  +  IYY»WD(  3,  4  )  ) /(HUU  )  *AREA) 
PR0D(1,5)=WD(1,4)*P1(4,5)+R(2)*WDC2,4)*P1(4,5) 
PR0D(1,6)=WD(1,3)*P1(3,6)+R(2)*WD(2,3)*P1(3,6) 

END  IF 

IF  (INDEX  .EQ.  21)  THEN 
CALL  C0EFP1(X,P1) 
CALL  VECTR(X,R) 

PR0D(1,1)=WD(2,1)»P1(2,1)+R(2)*WD(2,2)*P1(2,1) 
PR0D(2,1)=WD(3,1)»P1(3,2)+R(2)«WD(3,2)*P1(3,2) 
PR0D(3.1)=WD(4,1)*P1(4,3)+R(2)*WD(4,2)*P1(4.3) 
PROD(4,l)-(-IZZ*WD(3,4)  +  IYY*WD( 4, 3) )/ (HU( 1) *AREA) 
PR0D(5,1)=WD(4.1)*P1(4,5)+R(2)*WD(4,2)*P1(4.5) 
PR0D(6,1)»WD(3,1)»P1(3,6)+R(2)*WD(3,2)*P1(3,6) 

END  IF 

IF  (INDEX  .EQ.  22)  THEN 
CALL  C0EFP1(X,P1) 

PROD(l,l)=Pl(2,l)  *  WD(2,2)  *  Pl(2,l) 
PR0D(1,2)=P1(2,1)  *  WD(2,3)  *  Pl(3,2) 
PR0D(1,3)=P1(2,1)  *  WD(2,4)  *  Pl(4,3) 
PROD(l,4)=»0.D0 

PR0D(1,5)-P1(2,1)  *  WD(2,4)  *  Pl(4,5) 
PR0D(1,6)=P1(2,1)  *  WD(2,3)  *  Pl(3,6) 
PROD(2,l)=PROD(l,2) 

PR0D(2,2)=P1(3,2)  *  WD(3,3)  *  Pl(3,2) 
PR0D(2,3)=P1(3,2)  *  WD(3,4)  *  Pl(4,3) 
PROD(2,4)=0.D0 

PR0D(2,5)=P1(3,2)  *  WD(3,4)  *  Pl(4,5) 
PR0D(2,6)=P1(3,2)  *  WD(3,3)  *  Pl(3,6) 
PROD(3,l)=PROD(l,3) 
PROD(3,2)=PROD(2,3) 

PR0D(3,3)=P1(4,3)  *  WD(4,4)  *  Pl(4,3) 
PROD(3,4)«0.D0 

PR0D(3,5)=P1(4.3)  *  WD(4,4)  *  Pl(4,5) 
PR0D(3,6)=P1(4,3)  *  WD(4,3)  *  Pl(3,6) 
PROD(4,l)=0.D0 
PROD(4,2)=0.D0 
PROD(4,3)=0.D0 
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PROD(4,4)=(IZZ  *  WD(3,3)  +  IYY  *  WD(4, 3) ) / (MU( 1 ) "AREA) 

PROD(4,5)»0.D0 

PROD(4,6)=0.D0 

PROD(5,l)=PROD(l,5) 

PROD(5,2)=PROD(2,5) 

PROD(5,3)-PROD(3,5) 

PROD(5,4)=0.D0 

PR0D(5,5)=P1(4,5)  *  WD(4,4)  *  Pl(4,5) 

PR0D(5,6)=P1(4.5)  *  WD(4,3)  *  PH3.6) 

PROD(6,l)=PROD(l,6) 

PROD(6,2)*PROD(2,6) 

PROD(6,3)-PROD(3,6) 

PROD(6.4)=0.D0 

PROD(6,5)=PROD(5,6) 

PR0D(6,6)=P1(3,6)  *  WD(3,3)  *  PK3.6) 
END  IF 
IF  (INDEX  .EQ.  1)  THEN 

CALL  C0EFG1(X,G1) 

CALL  TRANSP(TERM,G1,4,6,G1T) 

CALL  STIFNS(C) 

CALL  MATMUL(PRD.C,4,4,G1,4,6,G1) 

CALL  HATMUL(PRD,G1T,6,4,G1,4,6,PR0D) 
END  IF 
IF  (INDEX  .EQ.  10)  THEN 

CALL  VECTR(X.R) 

CALL  TRANSP(TERM,R, 4,1. PROD) 
END  IF 

IF  (INDEX  .EQ.  20)  THEN 
CALL  COEFPKX.P1) 
CALL  TRANSP(TERH.P1.4.6,PROD) 
END  IF 
END 


SUBROUTINE  COEFG1  (X.G1) 

•  • 

*  THIS  SUBROUTINE  CALCULATES  THE  VALUES  FOR  THE  GAHA1  COEFFICIENT     * 


REAL»8  A ( 2 ) , D ( 2 ) , THETA ( 2 ) . ALPHA( 2 ) , AREA, E, G, 
+        IYY,IZZ,HU(2),GRAV, 
+        G1(4,6),L,X,K 

COHHON  /BLOCK1/  A,  D,  THETA, ALPHA, AREA, E.G. IYY, IZZ. MU.GRAV 

L=A(2) 
K  =  X/L 

DO  10  J=l,4 
DO  10  H=l,6 
10  G1(J,M)=0.D0 

Gl(l,l)  =  1/L 

Gl(2,3)  =  -(12. *K  +  6.)/L**2. 

Gl(2,5)  =  (6. *K  +  4.  )/L 

Gl(3,2)  =  GK2.3) 

GK3.6)  »  Gl(2,5) 

GK4.4)  *  GK1.1) 
END 
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SUBROUTINE    STIFNS    (S) 
************************************************************************ 

*  * 

*  THIS  SUBROUTINE  CREATES  THE  STIFFNESS  MATRICES  FOR  EACH  LINK     * 

*  GIVEN  THE  MATERIAL  AND  AREA  DATA  FROH  COMMON  BLOCK7  * 

*  * 
••A********************************************************************* 

REAL*8  A(2),D<2),THETA(2),ALPHA(2),AREA,E,G, 
+        IYY,IZZ,MU(2),GRAV,S(4,4),P 

COMMON  /BLOCK1/  A, D, THETA, ALPHA, AREA, E,G, IYY, IZZ,MU,GRAV 

DO  10  1=1,4 

DO  10  J=l,4 
10  S(I,J)=0.D0 

P  -  IYY+IZZ 

S(l,l)  »  E*AREA 

S(2,2)  -  E*IYY 

S(3,3)  =  E*IZZ 

S(4,4)  »  G*P 
END 


SUBROUTINE  VECTR  (X,R) 

*  • 

*  THIS  SUBROUTINE  CALCULATES  THE  POSITION  VECTOR  R  * 
»  * 

REAL*8  A(2),D(2),THETA(2),ALPHA(2),AREA,E,G, 
+        IYY,IZZ,MU(2),GRAV, 
+       R(4,1),L,X 

COMMON  /BLOCK1/  A, D, THETA, ALPHA, AREA, E,G, IYY, IZZ,MU,GRAV 

L=A(2) 

*  CREATING  THE  LOCAL  R  POSITION  VECTOR 

DO  10  J-1,4 
10      R(J,1)=0.D0 

R(l,l)  »  1.D0 
R(2,l)  =  X 

END 


SUBROUTINE  COEFP1  (X,P1) 

a************************************************************************ 

•  * 

*  THIS  SUBROUTINE  MILL  CALCULATE  THE  VALUES  FOR  THE  PHI1  COEFFICIENT  * 

«««*»«»***«««««•««*«*«««*«**««««*««»•••*«*•««»»»«*•«*«**■»*«**«***«•»«»**« 

REAL*8    A ( 2 ) , D ( 2 ) , THETA ( 2 ) , ALPHA ( 2 ) , AREA, E, G, 
+  IYY,IZZ,MU(2).GRAV, 

+  P1(4,6),L,X,K 

COMMON    /BLOCK1/    A, D, THETA, ALPHA, AREA, E,G, IYY, IZZ,MU,GRAV 

L=A(2) 
K    =   X/L 
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DO  10  J=l,4 

DO  10  M=1.6 
10  P1(J,H)=0.D0 

Pl(2,l)  ■   1.  +  K 
Pl(3,2)  »  -2.*K**3  -  3.*K»*2  ♦  1. 
Pl(4,3)  »   PK3.2) 
Pl(4,5)  =   L*K*(K**2  +  2."K  +  1.) 
Pl(3,6)  -   Pl(4,5) 
END 


SUBROUTINE  BINT  (TERM. XL. XU, Z, Y.NR.NC, INDEX, DD, HAT) 

*  GIVEN  THE  TEMPORARY  HATRICES  Z,  AND  Y  THE  INDEX  #  * 

*  CORRESPONDING  TO  INTEGRAL  DEFINITIONS,  THE  DIMENSIONS  OF  THE  * 

*  RESULTANT  MATRIX  PRODUCT  (ROM, COL),  AND  THE  W  FUNCTION  INDEX  (D)  * 

*  THIS  SUBROUTINE  WILL  INTEGRATE  THE  RESULTANT  PRODUCT  MATRIX  * 

*  OVER  THE  LENGTH  OF  THE  LINK  USING  GAUSSIAN  QUADRATURE  METHODS     * 

*  AN  ASSUMPTION  OF  CONSTANT  HEIGHT  DENSITY  AND  CROSS  SECTIONAL  * 

*  AREA  HAS  BEEN  MADE  WRT  THE  LINK/ACTUATOR  VOLUME  INTEGRAL.  * 

INTEGER  NR,NC, INDEX, TERM 

REAL*8  A(2),D(2),THETA(2),ALPHA(2),AREA,E,G, 
+        IYY,IZ2,HU(2),GRAV, 

+        AA,B,C,XL,XU,Z(NR,NC),Y(NR,NC).HAT(NR.NC),DD(4,4) 
COMMON  /BLOCK1/  A, D, THETA, ALPHA, AREA, E,G, IYY, I2Z,MU.GRAV 
AA=.5D0*(XU+XL) 
B=XU-XL 

C». 480 1449282487681 2D0»B 

CALL  INTPRD  (TERM, INDEX, DD, AA+C. NR, NC, Y) 
CALL  INTPRD  ( TERH, INDEX, DD,AA-C,NR, NC, Z) 
DO  10  1=1, NR 
DO  10  J=1,NC 
10  HAT(I,J)=.50614268145188130D-1*(Y(I,J)+Z(I,J)  ) 

C=».  398 33323870 68 1337D0*B 

CALL  INTPRD  (TERM, INDEX. DD, AA+C, NR, NC, Y) 
CALL  INTPRD  ( TERM, INDEX, DD,AA-C, NR,NC,Z) 
DO  20  1=1, NR 

DO  20  J=1,NC 
20  MAT(I,J)=MAT(I,J)+.11119051722668724D0*(Y(I,J)+Z(I,J)  ) 

C=. 26276620495816449D0»B 

CALL  INTPRD  (TERM, INDEX, DD, AA+C, NR,NC, Y) 
CALL  INTPRD  ( TERM, INDEX, DD, AA-C, NR, NC, Z ) 
DO  30  1=1, NR 

DO  30  J=1,NC 
30  MAT(I,J)=MAT(I,J)  +  . 15685 33 2293894 364D0* ( Y( I , J ) +Z ( I , J )  ) 

C=.917173212478  2490D-1*B 

CALL  INTPRD  ( TERM, INDEX, DD, AA+C, NR, NC, Y) 
CALL  INTPRD  ( TERM, INDEX, DD, AA-C, NR, NC, Z ) 
DO  40  1=1. NR 

DO  40  J=1.NC 
40        HAT(I,J)=BMMAT(I,J)  +  .18134189168918099D0*(Y(I,J)+Z(I,J)  )  ) 
DO  50  1=1, NR 

DO  50  J=1,NC 

IF(  INDEX.  EQ.  DTHEN 

HAT(I,J)=  MAT(I,J) 
ELSE 

MAT(I.J)=(MU(1)*AREA)  »  MAT(I,J) 
ENDIF 
50   CONTINUE 
RETURN 
END 
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SUBROUTINE  DFORM(H,K,UTERM,DC) 

*  * 

*  THIS  PROGRAM  CALCULATES  THE  DEFORHATION  TRANSFORMATION  MATRIX     * 

*  FOR  THE  KINETIC  AND  POTENTIAL  ENERGY  CONTRIBUTION  OF  POINT  LOAD   * 

*  • 

REAL*8  UTERM(6),DC(4,4) 

DO  10  1-1,4 
DO  10  J=l,4 
10  DC(I,J)=0.D0 

IF(M.EQ.1)THEN 

DC(1,1)»1.D0 

DC(2,2)=1.D0 

DC(3,3)=1.D0 

DC(4,4)=1.D0 
ENDIF 

IFUK.EQ.l)  .OR.  (K.EQ.7))THEN 

DC(2,1)=UTERM(1) 
ENDIF 

IF((K.EQ.2)  .OR.  (K.EQ.7) )THEN 

DC(3,1)=UTERH(2) 
ENDIF 

IF((K.EQ.3)  .OR.  (K.EQ.7) )THEN 

DC(4,1)-UTERM<3) 
ENDIF 

IF((K.EQ.4)  .OR.  (K.EQ.7))THEN 

DC(3,4)=- UTERM(4) 

DC(4,3)=UTERH(4) 
ENDIF 

IF((K.EQ.5)  .OR.  (K.EQ.7))THEN 
DC(2,4)=- UTERM(5) 
DC(4,2)=  UTERM(5) 
'  ENDIF 

IFUK.EQ.6)  .OR.  (K.  EQ.  7)  )THEN 

DC(2,3)=-UTERM(S) 

DC(3,2)=UTERM(6) 
ENDIF 
RETURN 
END 

SUBROUTINE  CC0EF( EPSLON, IOTA, H) 

A*********************************************************************** 

*  FORMULATING  THE  COEFFICIENTS  USED  IN  THE  SEQUENTIAL  INTEGRATION   * 

*  ADDITIONAL  DEFINITION  OF  VARIABLES  * 

*  Al-7 THE  COEFFICIENTS  FOR  NUMERICAL  INTEGRATION    * 

•ft********************************************************************** 

REAL*8  EPSLON, IOTA, H, 
♦        A0,A1,A2,A3,A4,A5,A6,A7 
COMMON  /BLOCK8/  AO, Al, A2, A3, A4, A5, A6, A7 
AO  »  l.D0/EPSLON/(H*H) 
Al  -IOTA/EPSLON/H 
A2  =  1.D0/EPSLON/H 
A3  =  0.5/EPSLON  -  1.D0 
A4  =  IOTA/ EPSLON  -  1.D0 
A5  =  H/2.D0  *  (A4  -  1 . DO ) 
A6  =  HM1.D0  -  IOTA) 
A7  =  IOTA  *  H 
END 


69 


SUBROUTINE  TRANSP  (  TERM, HI, R, CHIT) 

GIVEN  A  MATRIX  H  THIS  SUBROUTINE  WILL  RETURN  THE  TRANSPOSITION 
*      OF  THE  MATRIX,  MT  „ 

INTEGER  R,C,TERH 

REAL*8  M1(R,C),  M1T(C,R) 

DO  10  J  =  1,C 
DO  10  K=1,R 

M1T(J,K)=H1(K,J) 
10   CONTINUE 
END 

SUBROUTINE  RINT  ( I, H, NOME, J,K) 

»  * 

*  THIS  SUBROUTINE  WILL  PRINT  ANY  OF  THE  ARRAYS  USED,  IT  IS  * 

*  PRIMARILY  USED  FOR  CHECXING  THE  RETURNED  VALUES.  * 

*  * 

INTEGER  I,J,X,L,N,0 
REAL*8  M(I,J,X) 
CHARACTER*6  NAME,NOHE 

NAHE=»        • 

NAME=NOME 

DO  10  L  -  1,1 

WRITE(6,1000)  'ARRAY  ',NAME,'  FOR  I  ■  »,L 
DO  10  N  -  i,J 

WRITE(6,1100)  (H(L,N,0).0«1,X) 
10  CONTINUE 
1000  FORMAT  ( IX. A8, : AG, A9, II ) 
1100  FORMAT  (IX, 7(:  2X.E9.3)) 
*  END 


SUBROUTINE  REFORM 

*  REFORMULATING  THE  COEFFICIENT  MATRICES  OF  THE  SMALL  MOTION       * 

*  EQUATION  OF  HOTION  FOR  SOLUTION  BY  SEQUENTIAL  INTEGRATION         * 

*  • 

REAL* 8  MQQ, HQN ( 6 ) , MNN ( 6 , 6 ) , GN ( 6 , 6 ) , XN ( 6 , 6 ) , FN ( 6 ) , FQ, 
+        MN ( 6 , 6 ) , FNN ( 6 ) , FQTEMP, MQNTMP ( 6 )  , TEMP 1 (6,6), PRD (6,6) 

COHMON  /BLOCX5/  MQQ, HQN, MNN, GN, XN, FN, FQ 
COMMON  /BLOCX7/  MN.FNN 

FQTEMP=( 1 . DO /MQQ ) *FQ 
DO  10  1=1,6 
10      FNN(I)=FN(I)-MQN(I)«FQTEMP 
DO  20  1=1,6 

MQNTMP ( I ) = ( 1 . DO / MQQ ) *MQN ( I ) 
20   CONTINUE 

CALL  MATMUL(PRD,MQN,6,1, MQNTMP, 1,6, TEMPI) 
DO  30  1=1,6 

DO  30  J=l,6 
30  MN(I,J)=MNN(I,J)-TEMP1(I,J) 

END 
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* 


SUBROUTINE  TRANSF(TIME) 

*  THIS  SUBROUTINE  CREATES  THE  HARTENBERG-DENAVIT  MATRIX  FOR  EACH  * 
JOINT  AND  THE  TRANSFORMATION  MATRICES  RELATING  EACH  OF  THE  JOINT  * 
BODY  FIXED  COORDINATE  SYSTEMS  * 

*  ADDITIONAL  DEFINITION  OF  VARIABLES  * 
HD THE  HARTENBERG  DENAVIT  TRANSFORMATION  MATRIX  * 

REAL*8  A(2),D(2),THETA<2),ALPHA(2),AREA,E,G, 
+        I*TMZZ,MU(2),GRAV,HD(2,4.4),Q,QDOT,Q2DOT,TIME 
COMMON  /BLOCX1/  A,D,THETA, ALPHA, AREA, E,G, I YY, IZZ,MU,GRAV 
COMMON  /BLK15/  Q,QDOT,Q2DOT 
COMHON  /BLOCK2/  HD 

*  CREATING  THE  HARTENBERG-DENAVIT  MATRIX  FOR  EACH  I 

DO  10  1=1,4 

DO  10  J=l,4 

HD(1,I,J)=0.D0 
10  HD(2,I,J)=0.D0 

HD(1,1,1)  »   1.D0 

HD(1,2,2)  -   COS(Q) 

HD(1,2,4)  =■   SIN(Q) 

HD(1,3,2)  =   SIN(Q) 

HD(1,3,4)  =  -COS(Q) 

HD(1.4,3)  -   1.D0 

HD(2,1,1)  -   1. 

HD(2,2,1)  =   A(2)»COS(THETA(2) ) 

HD(2.2,2)  =   COS(THETA(2) ) 

HD(2,2,3)  =  -SIN(THETA(2)) 

HD(2,3,1)  »   A(2)*SIN(THETA(2) ) 

HD(2.3,2)  =   SIN(THETA(2) ) 

HD(2,3,3)  =   COS(THETA(2)) 

HD(2,4.4)  -   1.D0 
RETURN 
END 

SUBROUTINE  MATMUL  ( PRD, HAT1. R1.C1,HAT2, R2,C2, PROD) 

*  * 

*  THIS  PROGRAM  RECEIVES  TWO  MATRICES  TO  BE  MULTIPLIED  TOGETHER     * 

*  ENSURES  PROPER  DIMENSIONALITY  AND  RETURNS  THE  PRODUCT  MATRIX     * 

*  * 

INTEGER  R1,C1,R2,C2 

REAL*8  HAT1(R1,C1),MAT2(R2,C2),PR0D(R1,C2),PRD(R1.C2),SUM 

*  CHECKING  THE  DIMENSIONALITY  OF  THE  MATRICES  TO  BE  MULTIPLIED 

IF  (C1.NE.R2)  PRINT*, 'ERROR  IN  MULTIPLICTION  DIMENSIONS' 
»   MULTIPLYING  THE  MATRICES  TOGETHER 

DO  30  J=1,R1 

DO  20  K=1,C2 
SUM=0 . 0 
DO  10  M=1,C1 

IF(ABS(MAT1(J,M) ) .LE. l.E-30)  MATK J,H) =0 . DO 
CONTINUE 

IF (ABS(MAT2(M,K)).LE. l.E-30)  MAT2 (M, K ) =0 . DO 
CONTINUE 

SUM  =  SUM  +  MAT1(J,M)  »  HAT2(H,K) 
10  CONTINUE 

PRD(J,K)=SUM 
20      CONTINUE 
30   CONTINUE 

DO  50  1=1, Rl 
DO  50  J=1,C2 
50  PROD(I,J)=PRD(I,J) 

END 
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SUBROUTINE   WHATRX 

*  * 

*  THIS  SUBROUTINE  MULTIPLIES  THE  HARTENBERG-DENAVIT  MATRICES        * 

*  TOGETHER  TO  OBTAIN  THE  COORDINATE  TRANSFORMATION  MATRIX  * 

*  * 

*  ADDITIONAL  DEFINITION  OF  VARIABLES  * 

*  • 

*  HD THE  HARTENBERG  DENAVIT  TRANSFORMATION  MATRIX  » 

*  H TRANSFORMATION  MATRIX  * 

*  • 

REAL*  8  W(4,4),HD(2,4,4), a  SUM ( 4 , 4 ) , SUM 

COHMON  /BLOCX2/  HD 
COMMON  /BLOCX3/  W 

DO  10  L=l,4 

DO  10  H  =  l»4 

WSUM(L.H)=0.D0 
10  H(L,M)=0.D0 

DO  40  J  =  l,4 

DO  30  K  =  l,4 
SUM=0.0 
DO  20  H=»l,4 

SUM  -  SUM  +  HD(1,J,M)  *  HD(2,M,X) 
20  CONTINUE 

WSUM(J,X)=SUH 
30      CONTINUE 
40   CONTINUE 

DO  50  1=1,4 

DO  50  J=1.4 
50  W(I,J)=WSUH(I,J) 

END 
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APPENDIX  E 
FIGURES 
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Figure  1   N-ROSS  Terrestrial  Track 
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Figure  2   LFMR  Antenna  Subsystem  Configuration 
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Figure  3   The  Large  Motion  Position  e.  Constant  Torque, 
Zero  Tip  Mass. 
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Figure  4   The  Large  Motion  Position  9,  Constant  Torque, 
44  Kilogram  Tip  Mass. 
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Figure  5    The  Large  Motion  Velocity  9  Dot,  Constant 
Torque,  Zero  Tip  Mass. 
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Figure  6    The  Motion  Velocity  0  Dot,  Constant  Torque, 
44  Kilogram  Tip  Mass. 
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Figure  7   The  Large  Motion  Acceleration,  0  Double  Dot, 
Constant  Torque,  Zero  Tip  Mass. 
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Figure  8    The  Large  Motion  Acceleration  9  Double  Dot, 
Constant  Torque,  44  Kilogram  Tip  Mass. 
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Figure  9   The  Boom  Tip  Deflection  in  the  Axial  Direction, 
Constant  Torque,  44  Kilogram  Tip  Mass. 
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Figure  10  The  Boom  Tip  Deflection  in  the  Y  Transverse 
Direction,  Constant  Torque,  44  Kilogram  Tip 
Mass . 
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Figure  11  The  Boom  Tip  Deflection  in  the 
Direction,  Constant  Torque,  44 
Mass. 
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Figure  12 


The  Boom  Tip  Torsional  Deflection,  Constant 
Torque,  44  Kilogram  Tip  Mass. 
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Figure  13  The  Boom  Tip  Rotation  About 
Direction,  Constant  Torque, 
Mass. 
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Figure  14   The  Boom  Tip  Rotation  About  the  Local  Y 

Direction,  Constant  Torque,  44  Kilogram  Point 
Mass. 
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Figure  15  Simulation  Code  Flow  Chart 
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